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ABSTRACT 



This thesis consists of an interactive program that enables the student to study the 
orbital motion of satellites around the earth. The student can investigate the shape of 
a variety of orbits by varying the initial position and velocity of the satellite, or by sup- 
plying select orbital parameters i.e. initial orbital radius, eccentricity, and inclination. 
Satellite maneuvers can also be studied, like transfer orbits and inclination changes, by 
command velocity changes at any location in the orbit. Also the effects of the perturb- 
ing forces due to the oblateness of the earth, drag for low earth orbits, and gravitational 
attraction from the sun and moon can be investigated. The orbits are displayed in either 
the perifocal coordinate system around a model of the earth, or the ground track can 
be displayed on a map of the world. Orbital data is displayed below the orbital plot. 
The display is enabled by the use of display integrated software system and plotting 
language (DISSPLA) subroutines. 
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I. INTRODUCTION 



A visual aid for students new to orbital mechanics is required to comprehend fully 
the dynamics of orbital motion. This program is an interactive time step simulation 
program that calculates and plots cither unperturbed or perturbed elliptical orbits. The 
program interacts with the student in developing the initial orbit. Also the program 
enables the student with the ability to change the velocity of the satellite at a specific 
location in the orbit. This feature will permit the student to investigate the effects of 
commanded velocity changes as in perigee kicks, apogee kicks and inclination changes. 
The user can also modify the initial position and velocity of the satellite at the com- 
pletion of any orbit. 

The student is given an opportunity to investigate the elfccts of perturbing forces 
on the satellites orbit by choosing to have the program calculate the orbit with or with- 
out perturbing forces. The variation of parameters method, as seen in [Ref. 1: pp. 
396-407j. is used in calculating the perturbing orbit. The perturbing forces taken into 
consideration are the following: 

1. the oblateness of the earth 

2. drag for low earth orbits 

3. gravitational force of the moon 

4. gravitational force of the sun 

In order to review fully the operation of the program {included in appendix A ) and 
to uncover any problems or limitations that plagued the programming, the program has 
been divided up as follows: 

1. program design 

2. unperturbed orbit 

3. perturbed orbit 

4. velocity changes 

5. graphical plots 

The programming approach and equations used in each of the above sections will be 
examined in there respective chapters. A review of the coordinate systems used and their 



transformations between them are included in appendix B. Since all the equations used 
in the calculation of the orbital elements are from reference I. they will not be reviewed 
in each chapter but will be included in appendix C for a quick reference, equations from 
other sources will be referenced in their respective chapters. 

Examples of perturbed and unperturbed orbital plots for a variety of initial orbital 
parameters are included in appendix D. Included arc plots of low earth orbits, transfer 
orbits and geosynchronous orbits. 



II. PROGRAM DESIGN 



In designing this program an attempt was made to make it not only as user friendly 
as possible, but also to make the program as simple as possible to understand. To 
achieve these goals, the program would have to be written in a logical manner, in a 
computer language that is easy to follow, the program would have to run on terminals 
readily available to students (at the Naval Postgraduate School (NFS)), and the program 
would have to be easily used by students with a minimum amount of computer or orbital 
mechanics knowledge. 

FORTRAN was chosen as the programming language since it is a wildly used sci- 
entific language and it allows for very structured programming. By programming in a 
structured Format, the program can be expanded in the future with a minimum amount 
of time required to understand the programming code. FORTRAN also allows for 
double precision numbers to be used in the calculation of the orbit. This is critical when 
round off error in single precision could be greater then the actual change that one is 
trying to model. The equations in the descriptions of the program might not exactly 
match the equations in the listings because of special programming techniques which 
must be included in most computer programs to handle such problems as "division by 
zero". 

The display integrated software system and plotting language (DISSPLA) package 
available on the mainframe computer at NFS was used to enable a variety of graphical 
displays with a minimum amount of programming. DISSPLA has a set of subroutines 
that the programmer calls to display data contained in arrays. This requirement forces 
the program to load arrays with the satellites position in order for it to be plotted. The 
TEC 618 computer terminal and associative plotter was used for ease of gaining hard 
copy plots of the orbits and the diversity of locations that are available here at NPS. 
In order to run a program in DISSPLA the user must first define storage space of 1500k 
and designate temporary disk space, and then call DISSPLA with the program name. 
This is accomplished with the following commands: 

1. DEFINE STORAGE 1500K 
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2. I CMS 

3. TDISK 4 D1S 

4. DISSPLA ORBIT 



To make the program user friendly, the user is prompted for inputs via the keyboard. 
The entry is usually a number. A yes or no response can be entered by typing "V" or a 
"X". In most cases the program does a check to see if the input is appropriate. In order 
to make it as easy as possible for the student to get the desired orbit displayed, the 
program requires only the initial position and velocity of the satellite. The initial posi- 
tion and velocity of the satellite is supplied by the user in one of two ways. The user can 
input the position and velocity of the satellite, using the perifocal coordinate system 
(UK), or the user can lot the program place the satellite on the "l" axis of the UK system 
at the radius of perigee (RP) distance supplied by the user. This latter choice gives the 
initial location of the satellite, but to get the velocity the program will prompt the user 
for one of the following: 

1. the actual velocity in the UK system. 

2. the eccentricity (e) of the orbit. In which case the velocity is calculated from the 
following equations: 



a = 



IIP 

1 - e 



= semi-major axis 



EXR = - 



= energy mass 



Where u. = MG 

M = mass of earth 
G = Universal gravitational constant 



v = / 2( EXR 



u 

~RP^ 



3. the radius of apogee (RA). The velocity is calculated by first calculating the ec- 
centricity (e) from the following: 

_ RA - RP 
e RA -r RP 

With the eccentricity the same equations used above are used to calculate the ve- 
locity. 



In order to give the velocity a direction the inclination (i) of the orbit is required from 
the user. The following equations are used to calculate the velocity vector: 



v, = 0.0 



Vj - V cost i) 



’>> = v sini i) 

The program will check to ensure that the orbital eccentricity is less than 1.0, if it is not 
then the program will reject the inputs. After the initial input are accepted, the program 
will do calculations for the six orbital elements required to describe the size, shape and 
orientation of the orbit, and to pinpoint the position of the satellite along the orbit at a 
particular time. This classical set of six orbital elements are as follows: 

1. a. semi-major axis. 

2. e. eccentricity. 

3. i. inclination. 

4. £2. longitude of the ascending node. 

5. co. argument of perigee passage. 

6. T.time of perigee passage. 

The program actually calculates more orbital elements than the six classical elements 
required to plot the orbit, this is done in an effort to make the program as robust as 
possible. This will add in the ability to expand the program in the future. 

If the satellite is not initially at the perigee point then the satellite must first be 
stepped around to the perigee point. The program then enters a loop that calculates the 
orbit from the perigee point through one complete orbit around the earth and back to 
the perigee point. The orbit is calculated in steps of 2 times pi divided by an integer, i.e., 
2 times pi divided by 50. This step size was used to ensure a smooth orbit for display 
purposes and also to get within adequate distance to the perigee point or other location 
for a velocity change. After the loop is completed, the program will offer the user a 
choice of the following plots to check the orbit: 

1. perifocal 

2. groundtrack 

The program then goes into a loop offering the user the following choices until the user 
decides to end the program: 

1. plot another view of the same orbit. 
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III. UNPERTURBED ORBIT 



The subroutines that calculate the unperturbed orbit are the most widely used sub- 
routines in the entire program. These subroutines are called to step the satellite around 
to the perigee point from the user supplied initial position and velocity, to calculate the 
next unperturbed orbit, and for any velocity change. No matter which of these sources 
supply the initial position and velocity the program calculates the unperturbed orbit in 
the same manner. The only difference is where in the orbit the satellite is initially when 
these subroutines are called. Before the unperturbed subroutines are called, the orbital 
elements are calculated. 

The unperturbed subroutines are called by a single subroutine UNPRET' which has 
the following basic algorithm: 

1. Increment time by the time step size (DTE The time step was chosen as the period 
divided by fifty to give a smooth plot, but more importantly to ensure that the 
satellite is within an acceptable distance from a specific location for a velocity 
change. The angular error caused by the step size can be as much as PI 50 from 
the desired point for a circular orbit and will increase for more eccentric orbits. 
This error becomes a factor when the user is making velocity changes, and therefore 
it will be covered in that chapter in further detail. 

2. Calculate the new elements. The calculation of the new elements is the heart of this 
algorithm. The size, shape and orientation of the orbit remains unchanged. What 
is required is the position of the satellite along the orbit as a function of time. The 
problem becomes a matter to solve "the Kepler problem"-predicting the future po- 
sition and velocity of an orbiting object as a function of some known initial posi- 
tion and velocity and the time of flight [Ref. 1: p. 1S1J. An algorithm using these 
principles will follow: 

a. A time step (DT) is added to the time of flight(TF). time of flight is the elapsed 
time since the satellite passed the perigee point. 

TF = TF - DT 

b. The new mean anomaly (MA) is calculated from the new time of flight, and the 
mean motion (MM). 

MA = MM x TF 

c. With the new mean anomaly the new eccentric anomaly (EA) is calculated. 
Because the solution to the Kepler problem (MA = EA — ex sin(£Tj) is 
transcendental, an iterative solution based on the Newton method of root find- 
ing is used. The root in question is a solution to the equation 
(MA — EA — ex sinlE/f) = 0) . This algorithm takes the form of [Ref. 1: p. 222]: 

1 ) MA, = EA, — ex sin(£.£.) 



I : > ■ = E4 



Ml - MV\ 

( 1 — e x cos (EA n )) 



Where this equation is applied initially to EA 0 = MA and then reapplied 
until the difference between MA and MA„ becomes small enough to be ig- 
nored. 

d. The new true anomaly (v c ) is calculated from: 

cos -1 (<? - cos (EA )) 

'° ecos(£ff)— 1 

3. Calculate the new position and velocity. The position and velocity are calculated 
in the perifocal coordinate system (PQW). The PQW system uses the orbit as its 
fundamental plane and therefore requires only two coordinate to specify the satel- 
lite's position and velocity. The r„ coordinate is by definition always equal to zero. 
T he position of the satellite is calculated as: 

x w — r cos v 



Vl 



= 0 



The velocity of the satellite is calculated as: 

•\ = v -f ( “ sin v c^ 

= \'jr + cos v c ) 

v. = 0 

4. Store position and elements in arrays for plotting. In order for the program to plot 
the orbit the radius, true anomaly, inclination, and argument of perigee must be 
stored in arrays. The use of these arrays to plot the orbit will be explained in 
chapter 6. 

5. The process is repeated until the satellite is at the perigee point and the true 
anomaly is two pi. 

The procedure used to calculate the unperturbed orbit leave very little to be modified 
by a programmer. The only choices that had to be made concerned step size, how to tell 
the L’XPRCT subroutine that the perigee point had been reached, and a value of ac- 
ceptable error for newtons method. For the unperturbed orbit, the step size just had to 
be small enough to produce a smooth plot of the orbit. Two indicators for perigee were 
used, one was that the true anomaly was greater than 6.21 radians (two pi equals 6.2S 
radians) and the time from the previous perigee point will be greater then the period. 
The two indicators were logically 'and' together to ensure the perigee point was reached. 



The disparity between two pi and 0.21 radians is due to the error produced by the sat- 
ellite not beginning the orbit at exactly the perigee point and the step si/e used go 
around the orbit. The acceptable size of error for newtons method was set at l.<>T-lo. 
because for an unperturbed orbit this would be the major contributor to any error in the 
orbit and the magnitude of this error would be acceptable. However: in a perturbed 
orbit there are other factors contributing to determining the acceptable error, and these 
will be discussed in the next chapter. 
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IV. PERTURBED ORBIT 



The perturbed orbit uses the same basic routines as the unperturbed orbit in step- 
ping the satellite around the earth with one major difference, the perturbing forces 
produce a time rate of change of the orbital elements that must be applied at each time 
step. The variation of parameters method is used to determine this influence of the 
perturbing forces on the orbital elements. The analysis is simplified by using the orbital 
coordinate system RSW', as explained in appendix B. The basic algorithm is as follows 
[Ref. 1: p. 407]: 

1. At i — / 0 calculate six orbital elements. 

2. Compute the perturbing forces and transform it at t = r 0 to the 'RSW' SYSTEM. 

3. Compute the time rate-of-change of the elements. 

4. Calculate the change of elements for one time step, and add the changes to the old 
values at each step to get the new elements. 

5. From the new values of the orbital elements, calculate a position and velocity. 

6. Go to the step 2 and repeat until the final time is reached. 

The steps in the algorithm will be explained in the following sections: 

A. ORBITAL ELEMENTS 

The standard orbital elements a. e. i. £2, co and T (or M) will be used, where 
a = semi-major axis 
e = eccentricity 
i = inclination 

£2 = longitude of ascending node 
co = argument of perigee 
T = time of perigee passage 

(M c = mean anomaly at epoch = M — n(i — 0). The elements are calculated only 
at the beginning of the orbit from the initial position and velocity vectors. The elements 
are then changed continuously throughout the orbit by adding the changes due to the 
perturbing forces. For the perturbed orbit, the satellite will always begin at the perigee 
point. This is done so one complete orbit is from perigee point to perigee point. 
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B. COMPUTE PERTURBING FORCES 

The variation of parameters method requires that the perturbing forces be calculated at 
each step in the orbit. In order to do this a model of each perturbing force must be 
developed. The following perturbing forces where used in calculating the total perturb- 
ing force effecting the satellite: 

1. oblateness of the earth 

2. atmospheric drag 

3. gravitational attraction of the sun 

4. gravitational attraction of the moon 

The magnitudes of these forces have an enormous range of values and are dependent 
on the distance the satellite is from the perturbing body. Figure 1 on page 12 shows a 
graphical representation of the magnitude of the perturbing forces in a log-log plot of 
perturbing forces per unit mass [Ref. 2: p. IV-61]. The model of each of these forces 
follows: 

1. NON-SPHERICAL EARTH 

The earth is not perfectly spherical, but bulges around the equator. The polar 
and equatorial diameters are 12713.0 Km and 12756.3 Km. respectively. The oblateness 
results in a perturbing force per unit mass with these components in the 'RSW' coordi- 
nate system [Ref. 3: p. SI]: 

I —3 uJ-r~) 

F r = — (1-3 sin‘(zj sin'(R))) 






( -3u7u-;) 



( sin‘(z) sin(z/ 0 ) cos(> 0 )) 



Fu- 



( —3 uJ 2 i'j) 

r" 



( sin(z) cos(/) sin(zz 0 )) 



The variable and constants of these equations are defined below: 

1. Variables: 

a. zq = the argument of latitude and is equal to the true anomaly v 0 plus the ar- 
gument of perigee co. 

u , = v 0 - CJ 

b. r = the radius from the center of the earth to the satellite. 
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Figure 1. Comparison of perturbation magnitudes. 
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2. Constants: 

a. u = the gravitational parameter of the earth, 

(A/h 5 ) 

H = 398601.2 — r 1 - 

s~ 

b. J 2 = the second harmonic of oblateness coefficient, determined by experimental 
observations, 

J 2 = 1.0S23 E- 3 

c. r, = the mean radius of the earth, 
r e = 6.37S2£3A>?i 

2. ATMOSPHERIC DRAG 

The formulation of atmospheric drag equations are plagued with uncertainties 
of atmospheric fluctuations, frontal areas of orbiting object (if not constant), the drag 
coefficient, and other parameters. A fairly simple formulation will be given here. Drag, 
by definition, will be opposite to the velocity of the vehicle relative to the atmosphere. 
Thus, the perturbing force is 

F = — ( -vp— ) • CD • AR • DEX • vv 
2m 

The velocity vector is in the 'UK' system so the resulting force is also in the '1JK' sys- 
tem, Therefore a transformation to the 'RSW' system is required. 

The variables and constants of this equation are defined below: 

1. Variables: 

a. v = speed of vehicle. 

b. CD = the dimensionless drag coefficient. The drag coefficient CD has a value 
between 1 and 2. It takes a value near 1 when the mean free path of the at- 
mospheric molecules is small compared with the satellite size, and takes a value 
close to 2 when the mean free path is large compared with the size of the satel- 
lite. The drag coefficient will be modeled with CD = 2 when the satellites alti- 
tude is greater than 550km and equal to 1 otherwise. [Ref. 4: p. 295] 

c. DEN = atmospheric density at the vehicle's altitude. The density is spherically 
symmetric, and will be modeled using exponential steps using the parameters in 
Table 1 on page 14 and the following formula [Ref. 1: pp. 423-424]: 

(5(r) = . 
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Table 1. ATMOSPHERIC PAR AM ETERS AND YALl ES 



/ km) 


1 

5 | k 


z 


d(.-) 




1. 2 25 F -i)2 


4A4E-02 


0.0 


F2225E-02 








15o 


1.0E-03 


150-550 


1.79S4e>E-01 


4.3614E-02 


550 


3.0E-S 


550 

> 


1.0154S4E-07 


2.2169SE-07 


1500 


3.65E-09 








4100 


l.oE- 12 



2. Constants set to typical values: 

a. m = mass of the satellite, set equal to 100kg. 

h. AR = the cross-sectional area of the vehicle perpendicular to the direction of 
motion, set equal to 20m 2 

3. PERTURBING FORCE DUE TO HEAVENLY BODY 

The satellite will experience perturbation forces due to the gravitational effects 
of the sun and the moon. The perturbation force from a perturbing body is the differ- 
ence between the gravitational force due to the perturbing body at the satellite and the 
ravitational force the satellite would experience if it were at the center of the earth. 
From Figure 2 on page 15. the perturbing force per unit mass of the satellite is 

>' r L - ri r ti.i. 

f. = r — 

' t ^ . . \ 

\ r : , p - n r | ' ; 

The variable and constants are defined below: 

1 . Variables: 

a. r. = distance from the earth center for the perturbing body 

b. i p = unit vector from the earth to the perturbing body 

c. r = distance from earth center to the satellite 

d. /, = unit vector from the earth to the satellite 

2. Constants: 

a. jj.. = gravitational constant of the perturbing body = M p G 

The subscript p is to be replaced by s if the perturbing body is the sun. and by m if the 
perturbing body is the moon. We will assume that r < < r F then the equation above 
becomes 
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Fp = ( ~T )( 7~ )(3(/ r /p)/p - v) 

'/> p 

The unit vectors /, and i f can be written in terms of the 'UK' system as: 
i r = ( cos(Q) co$(«o) — sin(Q) cos(/) sin(tig))/ + ( cos(t^) sin(Q) + cos(w) cos(/) sin(t^)y + 

( sin(/) sin(uo))A' 

i p = ( cos(0,) cos(tijp) - sin(Q p ) cos(/ p ) sin(«o p ))7 + 

( cos(w 0p ) sin(Qp) + cos(co p ) cos (i p ) sin {u 0p ))J + ( sin(/ p ) sin(u 0 p))A' 

where £2. i, and Ug are the orbital elements of the satellites and fX., i p , and u ap are the 
orbital elements of the perturbing body. The formulas above use the 'IJK' system, and 
as such the resultant forces must be transformed to the 'RSW' system. Models of the 
sun and moon orbits are required to calculate r f and i p . The models used in the program 
for the sun and moon's orbits follows: (Ref. 3: pp. 73-74] 
a. S VIS’S POSITION 

In order to model the suns orbit, a number of simplifications had to be 
made in the actual parameters of the suns orbit. First the sun will be assumed to be in 
a circular orbit. This means that the radius (r) to the sun will be constant, and the ec- 
centricity (e) will equal 0.0 instead of its true value of 0.017. The other assumption will 
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be to place the sun on the T axis of the '1JK' system at the beginning of the program 
and have it progress through its orbit as the program runs. These changes will not effect 
the perturbing force in any noticeable magnitude. 

The following variables and constants where used in the program to model 
the suns orbit after applying the simplifications: [Ref. 3: pp. 75-78] 

1. Constants: 



a. Gravitational Constant: G = 6.6 7£ — 

b. Sun's Mass: m s = 1.99£30A'g 

c. Sun's Gravitational parameter: 



£ r = 






£20 



Mu- 



ll 



(Am : ) 

kg 2 



d. Sun's eccentricity: e. = 0.0 

e. Radius of orbit, assume sun is in circular orbit: r, = 1.49£ll/» 

f. Sun's inclination: si = 23.45 deg. = 4.09279709d-01 radians 

g. Longitude of ascending node: = 0.0 

h. Argument of perigee: cj s = 0.0 
2. Variables: 

a. The true anomaly of the sun's position as a function of the time the satellite has 
been in orbit: 



W 77 ) 356 x 24 x 3600 TT 

Where TT = true time, the time the satellite has been in orbit (sec) 

b. Sun's Position vector: r = r cos v 0; A 4- r sin v 0s Q 

c. Unit vector from the earth to the sun: i s = — 4— 



b. MOON'S POSITION 

In modeling the orbit of the moon, similar assumptions where used as with 
the sun. The moons orbit will be assumed to be circular, actually the eccentricity is 
equal to 0.055. By placing the moon initially on the T axis of the 'UK' system along 
with the sun, the gravitational forces of the two bodies will combine to a maximum. 
However; since the moons orbital period is only 27.3 days, the moon will not stay in this 
alignment and the magnitude of the combined forces will van.' with time. The inclination 
of the moons orbit is not constant, but drifts between IS. 3 and 28. 6 degrees in ten years. 
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Also the longitude of the ascending node (Q) oscillates between 13 and -13 degrees. To 
simplify this the inclination will be chosen as a constant 23.5 degrees and the longitude 
of the ascending node as 0.0 degrees. For the time period involved in calculating the 
perturbed orbit, these assumptions will not make any significant difference. 

The following variables and constants were used in the program to model 
the moons orbit, after applying the simplifications: 

1. Constants: 

(W) 

a. Gravitational Constant: G = 6.6/ E - 11 — — — 

kg' 

b. Moon's Mass: m„ = 7.2>5E22kg 

( Xnv) 

c. Moon's Gravitational Parameter: — GM m = 4.90E12 — 

kg 

d. Moon's eccentricity: e m = 0.0 

e. Radius of orbit, assume moon is in circular orbit: r m = 3.844£8£m 

f. Moon's inclination: i = 23.5deg. = 4. 10152374L- 1 radians 

g. Moon's longitude of ascending node: = 0.0 

h. Moon's argument of perigee: = 0.0 

i. Moon' 1 ' period: T = 27.3 days [period] 

2. Variables: 

a. The true anomaly of the moon's position as a function of the time the satellite 

has been in orbit: \\JTT) — — ■— — — : — TT 

^ ’ 2 .3 x 24 x 36UU 

b. Moon's position Vector: r = r cos v 0m P + r sin v 0m Q 

c. Unit vector from earth to moon: /_ = -4— 

I r m | 

The models of the sun and moons orbit calculates the position vector in the 'PQW' 
system and therefore the position vector must be transformed to the 'UK' system. 

C. RATE-OF-CHANGE OF ORBITAL ELEMENTS 

The derivations and equations of the rates-of-change of the orbital elements are 
contained in reference 1 pages 398 to 406. Therefore; only a summary of the actual an- 
alytic expressions for the rate-of-change of the parameters in terms of the perturbations 
will follow: 

1. Rate-of-change of the semi-major axis: 
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da _ r _ 
</7 L 



■]/•>[ 



Jc7 x 1 — e 



n r 



Vs 



Where n' is the mean motion of the satellites orbit. 



< / v 

" = \~ 
v a 

2. Rate-of-chance of the eccentricitv: 



de _ j- _\ 

dt _L 



1 — e~ sin v, 



n a 



•]/>+[ 



\ 1 - « 2 a 2 {\ - e 2 ) 



][' 



-Ws 



n a e 



. Rate-of-chanee of the inclination: 



at 



r cos 



y. 



\ 



4. Rate-of-change of the longitude of the ascending node: 
r sin u 0 



dt 



■y. 



n'a\ 1 — e" sin i 



5. Rate-of-change of the argument of perigee: 

({(•) , , , V) , 



Where. 



e cos v r 



( ) = r — 

1 dt ,r 1 n'ac 

(^l.-C^ltsin'od+T 



]/, 



C COS V, 



0 



■)y 



ii a ^ 1 — e 

6. Rate-of-chance of the eccentric anomaly; 



d e 



de 



dEA 



1 



[( sin v 0 + — )( 1 + e cos v 0 ) - ( cos v 0 + e)( — cos v 



dl 



di sin(E. l) [1 + e cos v c ]‘ 

7. Rate-of-chance of the mean anomalv: 



dMA dF. ! d, 



dt dl dl 



— sin! E. ! ) — e x cos 



(EA)dEA 

dl 



On 

dl 



(' - y 



IS 



e sin v 0 )] 



This equation reduces to the following for circular and ecliptic orbits 
(0 < = e < 1 ). 



dMA _ -1 r > 
di n'a L a 



(1 



cosv 0 ]/-;-[-^-][l 

n ae 



r n dn' 

: — ] sin 1 , /d — t — — 

n(l — e~ ) al 



Where the Rate-of-change of the mean motion: 
on' _ ~-M< , da 

ct l In' a* ci! 



[ref. 1 p. 396-407] 

D. NEW ORBITAL ELEMENTS 

The change of each element is calculated by multiplying the rate-of-change of the 
element by the time step ( DT). The change in the orbital elements arc then added to the 
current values of the elements to give the new orbital elements. With the new elements 
calculated, the satellite is stepped forward and the new position and velocity are calcu- 
lated in the same manner as the unperturbed orbit (chapter 3). Also as with the unper- 
turbed orbit, the process is repeated until the satellite is at the perigee point, indicated 
by the time of flight (TF) equal to the period of the perturbed orbit. 
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V. VELOCITY CHANGES 



The ability of the student to change the velocity of the satellite at any position in 
the orbit is a vital element in this program. With velocity changes the student can in- 
vestigate the effects of varying the satellites velocity as in transfer orbits and inclination 
changes. In order to simplify the program the unperturbed orbit is used throughout this 
routine. The velocity change algorithm used in the program follows: 

1. Rotate to velocity change location. 

The user is given the choice of changing the velocity of the satellite at the 
perigee, apogee or at any true anomaly. If the user chooses perigee or apogee as 
the change locations, the true anomaly is set equal to zero or pi radians respect- 
fully. With the location of the velocity change, the satellite is first stepped around 
to the desired true anomaly. The stepping is identical with the unperturbed orbit 
with the exception that the stepping terminates when the true anomaly is greater 
or equal to the desired true anomaly. With a step size of one fiftieth of the period, 
the satellite is actually stepped around to a location near the desired location. This 
variance can be reduced by decreasing the step size but this would increase the 
computation time. This error will be a major factor in precise calculations of 
transfer orbits, or any other orbital maneuver where precise velocity changes are 
required. However: this program is not a tool to calculate precise orbital maneu- 
vers. but rather a learning tool for the student to get a feel for the results of velocity 
changes in a satellite's orbit. 

2. Change the velocity. 

With the satellite at the desired location, the program calculates and displays 
for the user the satellite's current velocity, escape velocity and circular velocity (the 
velocity required to circularize the orbit). The program will not allow velocities 
greater than or equal to the escape velocity. The user is given the option to enter 
a new velocity in the UK' system or to change the magnitude of the velocity in the 
orbital plane. If the user chooses to change the velocity in the orbital plane, the 
program will prompt the user for the magnitude of the velocity change, and multi- 
ply this change by a unit vector in the direction of the satellites' velocity. This ve- 
locity change vector is then added to the satellites velocity vector, to calculate the 
new velocity vector. 

3. Calculate new elements. 

The orbital elements are calculated with the new velocity vector and the satel- 
lite's position vector. 

•4. Complete the orbit. 

The program will complete the orbit to the new perigee point using the satel- 
lite's position, new velocity and new elements. There are a number of problems 
that arise if the satellite is just stepped around to the perigee point. For example, 
with velocity changes in the orbital plane the apogee and perigee directions can 
physically swap. This is a problem when plotting with the perifocal coordinate 
s\ stem because the A', axis points toward perigee. To avoid problems like this the 
arrays used in plotting the orbit must be cleared and the satellite's current position 



and velocity be treated as initial conditions. However; to compare the old and new 
orbits there is a desire to retain as much of the previous orbit as possible. The 
velocity changes where divided into the following four cases to handle these prob- 
lems: 

a. Change velocity in the orbital plane at the perigee point with the new velocity 
greater than the circular velocity. The perigee point will remain the same so the 
satellite is stepped around using the unperturbed subroutines. 

b. Change velocity in the orbital plane at the perigee point with the new velocity 
less than or equal to the circular velocity. The perigee and apogee directions 
will switch so the plotting arrays are first cleared and stored with the current 
location data. Because the satellite is now at the apogee point the satellite is 
stepped around to the perigee point storing the second half of the orbit. The 
entire next orbit is calculated and stored to get a complete orbit. 

c. Change velocity in the orbital plane at the apogee point with the new velocity 
less than the circular velocity. The perigee and apogee directions will remain the 
same, so the satellite is stepped around to the perigee point completing the orbit. 

d. This last case catches all the following velocity changes; velocity change in the 
orbital plane at the apogee point with the new velocity greater than or equal to 
the circular velocity, velocity changes at any other true anomaly in the orbital 
plane, and any velocity change out of the orbital plane. T he plotting arrays are 
cleared and stored with the current location data. No matter where in the orbit 
the satellite is. the satellite is first stepped around to the perigee point, and to 
ensure a complete orbit is plotted the entire next orbit is also calculated and 
stored. 
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VI. GRAPHICAL PLOTS 



The program provides two types of graphical displays of the orbit, a display in the 
perifocal coordinate system and a display of the satellite's ground track. Each display 
type is useful in observing different aspects of the orbit. The perifocal display will allow 
the user to see how certain orbital parameters change with different initial positions and 
velocities, and also how the parameters change with velocity changes at varying posi- 
tions in the orbit. The ground track will enable the user to gain an appreciation for the 
physical location of the satellite above the earth, and see how the orbital parameter af- 
fects the path of the satellite. The ground track will also display the precession of a se- 
quence of orbits. Both displays plot the position steps to give the user an understanding 
of how the satellite speeds up at perigee and slows down around apogee. 

The DISSPLA package on the mainframe computer was used to enable the plotting 
of the orbits. The versatility of plotting subroutines of DISSPLA makes the actual 
programming of the orbit a simple matter of initializing DISSPLA for the type of mon- 
itor being used, setting up the plotting area, initializing the axis and axis scale, and then 
plotting the desired curve from points contained in arrays. This is a simplified explana- 
tion of DISSPLA, but for further details on DISSPLA programming refer to the 
DISSPLA user's manual [Ref. 5]. DISSPLA also supplies subroutines to draw a variety 
of projections of the world and fill the projections with coast lines, latitude lines and 
longitude lines. There are a couple of DISSPLA requirements that did require special 
handling in the program. The requirement that the data be supplied in arrays forced the 
program to load arrays with the required position and parameters and to keep a counter 
for the number in the arrays. The array format requires the size of the array be specified 
in the beginning of the program. The array size needs to be large enough to hold a 
number of orbits, but not so large as to waist storage space. The program will continue 
to add orbital data to the arrays until the user chooses to delete the previous orbits. If 
a new initial position and velocity is entered or if the arrays will overflow with the next 
orbit the arrays will automatically delete all previous orbits. DISSPLA also requires that 
all data be in single precision format. The program calculates all orbits in double pre- 
cision in order to limit the effect of round-off error, but by using the single precision data 
for plotting will not affect the accuracy of the plot in any way. 



The subroutines used to display the orbits will be covered in the following three 
sections; 

A. PERIFOCAL PLOT 

The plotting of the orbit in the perifocal coordinate system is the easier of the two 
types of plots. Since the perifocal coordinate system has the orbital plane as the fun- 
damental plane, the only requirements to describe the orbit in the perifocal coordinate 
system are arrays with the true anomaly and the radius to the satellite. To give the user 
a sense of the size of the plot, the axis length varies with the eccentricity and semi-major 
axis length. Also a plot of the earth is plotted to the same scale, with the pole or center 
of the plot on the origin of the axis. The latitude of the earth at the center of the plot 
will vary with the inclination of the orbit. This plot will allow the user to see a relative 
view of the satellite's coverage in the minus Z' axis direction of the perifocal coordinate 
system. 

B. GROUND TRACK 

The ground track plot is a very complex subroutine compared with the perifocal 
plot. Because the ground track is not a continuous curve a procedure to handle the 
satellite ending at one end of the plot and wrapping around to the other end was devel- 
oped. The wrap around problem is avoided in most orbits by plotting the orbit in seg- 
ments with the following two rules. Each segment begins at the beginning of a new plot 
or at the edge of the plot area, and ending when the satellite would wrap around to the 
other side of the plot. At the beginning of a segment if the position of the satellite is 
within five degrees of the edge of the plot, that position and any other positions within 
that five degree boundary will not be plotted. The segment will end when the satellite 
is within ten degrees of the edge of the plot. The above restrictions imposed on the 
segments of the plot will not substantially affect the interpretation or usefulness of the 
plot. The ground track is plotted on top of a cylindrical equidistant projection of the 
world, with the world coast lines and a longitude-latitude grid for reference. 

C. DATA 

Information concerning the orbit is displayed on the lower half of the plot. The 
information is designed to supply the user with enough of the basic orbital elements and 
other parameters affecting the orbit to be able to evaluate what basic type of orbit the 
satellite is in. and the effects of velocity changes and perturbing forces have on the orbit. 
The following data are plotted: inclination! i >. semi-major axis (a), eccentricity (e). period 
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(per), c 
emerus 



ipc r ee and perigee velocity and radius, average time rate-of-chanee of orbital el- 
. and tae average magnitude of perturbing forces per unit mass. 



MI. CONCLUSIONS AND RECOMMENDATIONS 



The program supplies the student with an interactive tool to study the orbital mo- 
tion of satellites around the earth. The student can investigate a variety of orbits by 
varying the orbital parameters, command velocity changes, and observe the effects of 
perturbing forces. 

The student is provided with two options for entering the initial position and veloc- 
ity of the satellite. The program could be expanded to provide the student with the ad- 
ditional options of entering either orbital parameters or a ground observation data and 
have the program calculate the initial position and velocity from this data. Also the 
student is limited to orbits with eccentricities less than one (elliptic orbits). The program 
could be also be expanded to include more eccentric orbit for Lunar, interplanetary, and 
missile trajectories. The perturbing orbit is calculated for orbits around the earth with 
relatively small perturbing forces in relation to the earths gravitational force. This fact 
will cause the program to produce false results if the student tries to calculate lunar 
trajectories. Special routines would have to be employed when the perturbing force (the 
moons gravitational attraction) is comparable to the earths gravitational attraction. 
This will not become a factor for studying current satellite orbits out to the 
geosynchronous radius of 42241.1km. 

The velocity change subroutines move the satellite to a location close to the desired 
location before a velocity change is imposed. By reducing the step size in the velocity 
change subroutine, this error could be reduced. Precise orbital transfer maneuvers can 
be modeled by reducing this error caused by the positioning of the satellite prior to 
changing the velocity. The program will currently provide the student with useful plots 
for gaining experience with various transfer orbits by varying the magnitude and location 
of the velocity changes. 

The output of the calculations of the orbit are arrays loaded with the satellite's po- 
sition and select orbital parameters. The DISSPLA subroutines that plot the points are 
not unique. The program would become portable to personal computers with these 
graphics subroutines written in FORTRAN and included in the program. 

A final recommendation is that the display of the ground track could be modified 
to show ground coverage, number of satellites in a constellation, and other elements 
necessary for planning a real-world artificial satellite application. 



APPENDIX A. ORBIT PROGRAM 

PROGRAM ORBIT ORBOOOIO 

THIS PROGRAM IS AN INTERACTIVE TIME STEP SIMULATION OF ORB00020 

SATELLITES AROUND THE EARTH. PERTURBED AND UNPERTURBED ORBITS ORB00030 

ARE CALCULATED AND PLOTTED. VELOCITY CHANGES ARE ALSO PERMITTED ORB00040 



AT SPECIFIED TRUE ANOMALIES. 




ORB00050 










ORB00060 


A LIST 


OF VARIABLES USED BY THE MAIN PROGRAM 


FOLLOWS: 


ORB00070 


A 


= 


SEMI -MAJOR AXIS 




ORB00080 


AL 


= 


ARGUMENT OF LONGITUDE 




ORB00090 


AP 


= 


ARGUMENT OF PERIGEE 




ORBOOIOO 


CHTA 


= 


VELOCITY CHANGE LOCATION TRUE ANOMALY 




ORBOOllO 


DT 


= 


TIME STEP 




ORB00120 


E 


= 


ECCENTRICITY 




ORB00130 


EA 


= 


ECCENTRIC ANOMALY 




ORB00140 


El 


= 


I VECTOR OF ECCENTRICITY 




ORB00150 


EJ 


= 


J VECTOR OF ECCENTRICITY 




ORB00160 


EK 


= 


K VECTOR OF ECCENTRICITY 




ORB00170 


FR 


= 


R VECTOR OF TOTAL FORCE 




ORB00180 


FS 


= 


S VECTOR OF TOTAL FORCE 




ORB00190 


FW 


= 


V VECTOR OF TOTAL FORCE 




0RB00200 


H 


= 


ANGULAR MOMENTUM 




ORB00210 


HI 


= 


I VECTOR OF ANGULAR MOMENTUM 




ORB00220 


HJ 


= 


J VECTOR OF ANGULAR MOMENTUM 




ORB00230 


HK 


= 


K VECTOR OF ANGULAR MOMENTUM 




ORB00240 


I 


= 


INCLINATION 




ORB00250 


IOPTl= 


PERTURBED OR UNPERTURBED OPTION 




ORB00260 


i opt: 


i — 


OPTIONS: PLOT NEXT ORBIT. CHANGE INITIAL VALUES. 


ORB00270 


CHANGE 


VELOCITY, PLOT ANOTHER VIEW OF ORBIT, 


QUIT 


ORB00280 


LAN 


= 


LONGITUDE OF ASCENDING NODE 




ORB00290 


LP 




LONGITUDE OF PERIGEE 




ORB00300 


MA 


= 


MEAN ANOMALY 




ORB00310 


MM 


= 


MEAN MOTION 




ORB00320 


MU 


= 


GRAVITATIONAL PARAMETER 




ORB00330 


N 


= 


ASCENDING NODE 




ORB00340 


NT 


= 


I VECTOR OF ASCENDING NODE 




ORB00350 


NJ 


= 


J VECTOR OF ASCENDING NODE 




ORB00360 


NK 


= 


K VECTOR OF ASCENDING NODE 




ORB00370 


NUM 


= 


STEP COUNTER 




ORB00380 


P 


= 


SEMI-LATUS RECTUM 




ORB00390 


PER 


= 


PERIOD OF ORBIT 




ORB00400 


PI 


= 


PI 




0RB00410 


RA 


= 


RADIUS OF APOGEE 




ORB00420 


RE 


= 


RADIUS OF EARTH 




ORB00430 


R 


= 


ORBITAL RADIUS 




0RB00440 


RI 


= 


I VECTOR OF ORBITAL RADIUS 




0RB00450 


RJ 


= 


J VECTOR OF ORBITAL RADIUS 




0RB00460 


RK 


= 


K VECTOR OF ORBITAL RADIUS 




0RB00470 


T 


= 


TIME COUNTER IN ORBIT 




0RB00480 


TA 


— 


TRUE ANOMALY 




ORB00490 


TDA 


= 


TOTAL CHANGE IN SEMI -MAJOR AXIS 




ORB00500 


TDAP 


= 


TOTAL CHANGE IN ARGUMENT OF PERIGEE 




ORB005 10 
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TDE = TOTAL CHANGE IN ECCENTRICITY 

TDK = TOTAL CHANGE IN ANGULAR MOMENTUM 

TDI = TOTAL CHANGE IN INCLINATION 

TDMA = TOTAL CHANGE IN MEAN ANOMALY 

TDMM = TOTAL CHANGE IN MEAN MOTION 

TDLAN= TOTAL CHANGE IN LONGITUDE OF ASCENDING NODE 

TF = TIME OF FLIGHT 

TFDRA= TOTAL FORCE OF DRAG 

TFEA = TOTAL FORCE OF EARTH’S OBLATENESS 

TFMO = TOTAL FORCE FROM MOON 

TFSU = TOTAL FORCE FROM SUN 

TL = TRUE Longitude AT EPOCH 

TT = TRUE TIME SINCE SATELLITE HAS BEEN IN ORBIT 

V = SATELLITE VELOCITY 

VI = I VECTOR OF SATELLITE VELOCITY 

VJ = J VECTOR OF SATELLITE VELOCITY 

VK = K VECTOR OF SATELLITE VELOCITY 

A LIST OF THE ARRAYS USED FOLLOWS: 

AINRAY = INCLINATION 
APRAY = ARGUMENT OF PERIGEE 
RARAY = RADIUS 
RIRAY = I VECTOR OF RADIUS 
RJRAY = J VECTOR OF RADIUS 
RKRAY = K VECTOR OF RADIUS 
TARAY = TRUE ANOMALY 
TIMRAY = TIME 

A LIST OF SUBROUTINES CALLED BY THE MAIN PROGRAM WILL FOLLOW: 
CALCEL = CALCULATES THE ORBITAL ELEMENTS 

CHGVEL = ALLOW THE USER TO CHANGE THE VELOCITY OF THE SATELLITE 

INPUTS = PROMPTS USER FOR INITIAL POSITION AND VELOCITY 

INTSUM = INITIALIZES THE SUMS IN THE ARRAYS 

NEWELT = CALCULATE NEW ORBITAL ELEMENTS FROM TIME STEP 

NEWPOS = CALCULATE NEW POSITION VECTOR 

NEWVEL = CALCULATE NEW VELOCITY VECTOR 

OPTION = GIVE THE USER THE OPTIONS Permitted IN THE PROGRAM 

PLOTS = PLOTS THE ORBITS 

PRETUR = CALCULATES THE PERTURBED ORBIT 

STORE = STORE THE POSITION DATA IN ARRAYS 

UNPRET = CALCULATE THE UNPERTURBED ORBIT 

BEGIN MAIN PROGRAM 

DOUBLE PRECISION PI ,MU,RI ,RJ,RK,R,VI ,VJ,VK,V,HI ,HJ,HK,H, 

+ NI ,NJ,NK ,N ,P ,EI ,EJ,EK,E , A , I , LAN , AP ,TA , AL, LP ,TL, PER ,EA , 

+ MM , MA , T , DT , TF , FR ,FS , FW ,TT , CHTA , RA , VA , TEMPTA , RE 

DIMENSION TARAY(500) , RARAY (500) ,RIRAY(500) ,RJRAY(500) ,RKRAY(500) , 
+ AINRAY(500) ,APRAY(500) ,TIMRAY(500) 

CHARACTER* 1 , LOOP , YORN , ORLOOP 

PI = 3. 141592653589794 



0RB00520 
ORB00530 
ORB00540 
ORB00550 
0RB00560 
ORB00570 
ORB00580 
ORB00590 
ORB00600 
ORB006 10 
ORB00620 
ORB00630 
0RB00640 
0RB00650 
ORB00660 
ORB00670 
ORB00680 
0RB00690 
ORB00700 
ORB007 10 
ORB00720 
ORB00730 
ORB00740 
ORB00750 
ORB00760 
ORB00770 
0RB00780 
0RB00790 
ORB00800 
ORB00810 
ORB00820 
ORB00830 
0RB00840 
ORB00850 
ORB00860 
ORB00870 
0RB00880 
ORB00890 
ORB00900 
ORB009 10 
0RB00920 
ORB00930 
ORB00940 
ORB00950 
ORB00960 
ORB00970 
ORB00980 
ORB00990 
ORBOIOOO 
ORBOIOIO 
ORBO 1020 
ORB01030 
ORB01040 
ORBO 1050 
0RB01060 
ORB01070 



MU = 3. 986012D+05 
RE = 6. 378145D+03 

USER INTRO TO PROGRAM 
CALL INTRO 

ENTERED MAIM PROGRAM LOOP 
LOOP = 'Y' 

10 IF (LOOP . EQ. ’ Y ' ) THEM 

INITIALIZE STEP COUNTER AMD TRUE TIME 
20 NUM = 1 

TT = 0. 0 

PROMPT USER FOR INITIAL POSITION AND VELOCITY 
CALL INPUTS(RI ,RJ ,RK,R, VI , VJ, VK, V,MU ,LOOP ,PI ) 

EXIT PROGRAM 
IF (LOOP .EQ. 'N') THEN 
GOTO 10 

END IF 

CALCULATE AND STORE ORBITAL ELEMENTS 

CALL CALCEL(RI , RJ , RK,R, VI , VJ, VK, V,EI ,EJ,EK,E , A, I , LAM, 

+ LP ,TA , PER,EA ,MA , AP , AL,TF , P, PI ,MU,MM,N,H,HI ,HJ) 

CALL STORE (RI ,RJ ,RK ,R ,TA , RIRAY , RJRAY , RKRAY , RARAY , TARAY , 

+ NUM , I , AP , A INRAY , APRAY , TT , TIMRAY) 

PRINT DATE FOR USER TO REVIEW 
PRINT*, 'VI =' , VI, ' KM/S' 

PRINT*, 'VJ =' , VJ, ' KM/S’ 

PRINT*, 'VK =' , VK, ' KM/ S' 

PRINT*, ' V =' , V, ' KM/ S' 

PRINT*, 'RI =' , RI,' KM' 

PRINT*, 'RJ =' , RJ, ' KM' 

PRINT*, 'RK =' , RK, ' KM' 

PRINT*, ' R =' , R, ' KM' 

PRINT* , ' ECCENTRIC ITY = ' , E 
DEGI = SNGL(( 180. 0/PI )*I) 

PRINT*, ’ INCLINATION =' ,DEGI , ' DEGREES ' 

PERHRS = SNGL( PER/3600. 0) 

PRINT* , ' PERIOD =' , PERHRS , ’ HOURS ' 

PRINT*, 'ARE THESE VALUES CORRECT? ENTER "Y" OR "N" :' 
READ*, YORM 

CALL EXCMS( ' CLRSCRN' ) 

IF (.NOT. YORN . EQ. ' Y ' ) THEN 
GOTO 20 

END IF 

CALCULATE TIME STEP AND SET TIMER TO ONE TIME STEP 
DT = PER/ 50 
T = DT 

STEP SATELLITE TO PERIGEE POINT AND RECORD 
50 IF ((TA. GT. 0.063). AND. (TA. LT. 6. 21)) THEN 

TT = TT + DT 



ORB01080 
ORB01090 
0RB01100 
ORBOlllO 
ORB01120 
ORB01130 
ORB01140 
ORB01150 
ORB01160 
ORB01170 
ORB01180 
ORB01190 
ORB01200 
ORB01210 
ORBO 1220 
ORB01230 
ORBO 1240 
ORB01250 
ORB01260 
ORB01270 
ORBO 1280 
ORB0129O 
0RB01300 
ORB01310 
ORB01320 
ORB01330 
ORB01340 
ORB01350 
ORB01360 
ORB01370 
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CALL XEVELK MM , MA , E , E A , TA , TF , DT , P I , PER ) 

CALL NP0S( RI , RJ , RK ,R, LAN , AP , I , TA.A.E) 

CALL NVEL( E , P , TA , LAN’ , AP , I , VI , VJ , VK , V , MU ) 

NUM = SUM + 1 

CALL STORE ( RI , RJ , RK , R , TA , RIRAY , RJRAY , RKRAY , RARAY , TARAY , 
NUM , I , AP , A INRAY , APRAY , TT , TIMRAY ) 

T = T + DT 
GOTO 50 

END IF 

CALCULATE ELEMENTS FROM PERIGEE POINT 

CALL CALCEL(RI ,RJ,RK,R,VI ,VJ,VK,V,EI ,EJ,EK,E,A, I , LAN, 

LP ,TA,PER,EA,MA, AP , AL,TF , P ,PI ,MU ,MM ,N,H,HI ,HJ) 

DT = PER/ 50 
T = DT 

STORE FIRST Unperturbed ORBIT 

CALL UNPRET(DT , PER , AL , LAN , AP , I ,RI ,RJ ,RK ,R, VI , VJ, VK, V, 

MU , PI ,H , A, E ,N,TA,P , MM, MA,EA,TF ,T, NUM , RIRAY, RJRAY, 
RKRAY , RARAY , TARAY , A I NRAY , APRAY , T I MRAY , TT ) 



0RB01870 
0RB01880 
0RB01890 
0RB01900 
0RB019 10 



BEGIN NEW ORBIT OPTIONS 
I0FT1 = 1. Unperturbed ORBIT 
= 2. Perturbed ORBIT 
= 3. QUIT 

I0PT2 = 1. PLOT NEXT ORBIT 

= 2. CHANGE INITIAL VALUES 

= 3. CHANGE VELOCITY AT A SPECIFIC TRUE Anomaly 
= 4. PLOT ANOTHER VIEW OF SAME ORBIT 

ALSO ASKED IF WANT TO CLEAR ALL PREVIOUS ORBITS 

CALCULATE ELEMENTS AT PERIGEE 

CALL CALCEL(RI , RJ,RK,R,VI , VJ, VK,V,EI ,EJ,EK,E, A, I , LAN, 

LP,TA,PER,EA,MA,AP,AL,TF,P,PI,MU,MM,N,H,HI,HJ) 

CHECK FOR POSSIBLE ARRAY OVERFLOW 
IF (NUM . GT. 425) THEN 

PRINT*, 'ARRAYS ARE FULL' 

PRINT*, 'PREVIOUS ORBITS WILL BE ERASED! ' 

NUM = 1 

CALL STORE ( RI , RJ , RK , R , TA , R IRAY , RJRAY , RKRAY , RARAY , TARAY , 
NUM , I , AP , A I NRAY , APRAY , TT , TIMRAY) 

ENDIF 

PROMPT USER FOR DESIRED OPTION 
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CALL OPT I ON' ( IOPT1 , IOPT2 , NUM , RIRAY ,RJRAY , RKRAY , RARAY , 

TARAY , A INRAY , APRAY .TIMRAY) 

Initialize SUMS FOR FORCE AND ORBITAL ELEMENT CHANGES TO ZERO 
CALL INTSUMC TFEA , TFSU , TFMO , TFDRA , TD I ,TDA , TDE , TDMM , TDMA , TDLAN , 
TDH, TDAP) 

SET TIME COUNTER TO ONE TIME STEP 
T = DT 

OPTION: PLOT THE NEXT ORBIT 
IF ( IOPT2 .EQ. 1) THEN 

CALCULATE AND PLOT UNPERTURBED ORBIT 
I F C IOPT1 .EQ. 1) THEN 

CALL UNPRET( DT , PER , AL , LAN , AP , I , RI , R J , 

RK , R , V I , VJ , VK , V , MU , P I , H , A , 

E, N, TA , P, MM, MA,EA,TF,T, NUM, RIRAY, RJRAY, RKRAY, 
RARAY , TARAY , A I NRA Y , APRAY , T IMRAY , IT ) 

CALL PLOTS (RIRAY , RJRAY , RKRAY , RARAY , TARAY , NUM , 
PI,I,LP,A,E,TF , AINRAY , APRAY , TIMRAY , 

TFEA , TFSU , TFMO , TFDRA , PER , 

TDI , TDA , TDE , TDMM , TDMA , TDLAN , TDH , TD AP , 

MM , MA , LAN , H , AP , R , V ) 

CALCULATE AND PLOT PERTURBED ORBIT 
ELSEIF( IOPT1 .EQ. 2) THEN 

CALL PRETURC DT , PER , AL , LAN , AP , I , 

RI,RJ,RK,R,VI,VJ,VK,V,FR,FS,FW, 

MU , P I , H , A , E , N , TA , P , MM , MA , E A , TF , T , NUM , 

RIRAY , RJRAY , RKRAY , RARAY , TARAY , AINRAY , APRAY , 
TIMRAY , TT , TFEA , TFSU , TFMO , TFDRA , 

TDI , TDA , TDE , TDMM , TDMA , TDLAN , TDH , TDAP ) 

CALL PLOTS (RIRAY , RJRAY , RKRAY , RARAY , TARAY , NUM , 

P I , I , LP , A , E , TF , AINRAY , APRAY , TIMRAY , 

TFEA , TFSU , TFMO , TFDRA , PER , 

TD I , TDA , TDE , TDMM , TDMA , TDLAN , TDH , TDAP , 

MM , MA , LAN , H , AP , R , V) 

END IF 

GOTO THE BEGINNING OF THE PROGRAM TO CHANGE THE INITIAL VALUES 
ELSEIF ( IOPT2 . EQ. 2) THEN 
GOTO 20 

CHANGE VELOCITY AT A SPECIFIC TRUE ANOMALY AND 

PLOT THE NEW ORBIT 

ELSEIF ( IOPT2 . EQ. 3) THEN 

CALL CHGVEL(DT, PER, AL, LAN, AP ,I,RI,RJ,RK,R, 

VI ,VJ, VK,V,MU,PI , 

H , A , E , N , TA , P , MM , MA , EA , TF , T , NUM , RIRAY , 
RJRAY , RKRAY , RARAY , TARAY , AINRAY , APRAY , 
TIMRAY, TT, El, EJ,EK,LP, HI, HJ.IOPT1, 

TFEA , TFSU , TFMO , TFDRA , TD I , TDA , TDE , TDMM , 
TDMA , TDLAN , TDH , TDAP ) 

CALL PLOTS ( RIRAY , RJRAY , RKRAY , RARAY , TARAY , NUM , 

PI , I , LP , A , E , TF , AINRAY , APRAY , TIMRAY , 
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+ TFEA , TFSU , TFMO , TFDRA , PER , 

+ TDI , TDA , TDE , TDMM , TDMA , TDLAN , TDH , TDAP , 

4- MM , MA , LAN , H , AP , R , V ) 

PLOT ANOTHER VIEW OF THE SAME ORBIT 
ELSEIF ( IOPT2 . EQ. 4) THEN 

CALL PLOTS( RIRAY , RJRAY , RKRAY , RARAY ,TARAY , NUM , 

+ PI , I , LP,A,E,TF,AINRAY,APRAY,TIMRAY, 

+ TFEA, TFSU, TFMO, TFDRA, PER, 

+ TD I , TDA , TDE , TDMM , TDMA , TDLAN , TDH , TDAP , 

+ MM , MA , LAN , H , AP ) 

* STOP THE PROGRAM 

ELSEIF ( I0PT2 . EQ. 5) THEN 
GOTO 90 

ELSE 

PRINT*, 1 INVALID ENTRY!' 

GOTO 80 

END IF 

* CHECK IF SATELLITE Impacted THE EARTH AND GO TO THE BEGINNING 
IF (R . LE. 6450.0) THEN 

PRINT*, 1 SATELLITE WILL IMPACT THE EARTH!!! ' 

PRINT*, 'PROGRAM WILL RESET TO THE BEGINNING' 

GOTO 20 

END IF 

* GOTO THE TOP OF THE OPTION LOOP 
GOTO 80 

* GIVE THE USER A CHANCE TO RECOVER THE PROGRAM 

90 PRINT*, 'THIS IS YOUR LAST CHANCE! ' 

PRINT*, 'DO YOU WANT TO CONTINUE?' 

PRINT*, 'AND GOTO THE Beginning OF THE PROGRAM?' 

PRINT*, 'ENTER "Y" OR "N" : ' 

READ* , LOOP 
PRINT*, LOOP 
GOTO 10 
END IF 

* DISSPLA SUBROUTINE TO TELL GRAPHICS TERMINAL PLOTTING 
SESSION IS DONE 

CALL DONEPL 

STOP 

END 

***************************************** 



SUBROUTINE INTRO 

THIS SUBROUTINE WILL GIVE THE USER A Brief INTRODUCTION OF THE 
USES OF THE PROGRAM 

PRINT*, ' THIS PROGRAM IS A GRAPHICS DISPLAY OF Satellite ORBITS. ' 
PRINT*, 'YOU WILL BE ASKED TO INPUT THE INITIAL VELOCITY AND' 
PRINT*, 'POSITION VECTORS OF THE Satellite. THE PROGRAM WILL ' 
PRINT*, 'THEN CALCULATE THE ORBITAL PARAMETERS AND THE ' 
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PRINT*, 

PRINT*, 

PRINT*, 

PRINT*, 

PRINT*, 

PRINT* , 

PRINT*, 

PRINT* , 

PRINT* , 

PRINT*, 

PRINT*, 

PRINT*, 

PRINT*, 

PRINT*, 

PRINT*, 

PRINT*, 

PRINT*, 

PRINT*, 

PRINT*, 

RETURN 

END 



Unperturbed ORBIT. THE USER WILL THEN HAVE THE' 

CHOICE OF DISPLAYS: ' 

-PERIFOCAL (SHOWS RELATIVE SIZE OF ORBIT)’ 
-Equatorial (SHOWS ORBIT INCLINED, USER INPUT’ 
LONGITUDE TO VIEW AT)' 

-GROUND TRACK' 

I 

THE USER IS THEN ASKED TO CHOOSE ONE OF THE FOLLOWING: ' 
-Unperturbed ORBITS' 

-Perturbed ORBITS’ 

-VELOCITY CHANGES' 

THE USER"S CHOICE WILL BE USED IN DEVELOPING THE' 
GRAPHICAL OUTPUT. ’ 

T 

THE USER IS THEN GIVEN THE FOLLOWING CHOICES: ’ 

-CLEAR ALL THE PREVIOUS ORBITS' 

-CHANGE THE INITIAL PARAMETERS' 

-CHANGE VELOCITY AT A SPECIFIC TRUE Anomaly' 

-PLOT ANOTHER VIEW OF THE SAME ORBIT' 



* * * -v * * -A- * -,V * *********** * * ****** * * * * * ** * -,V ** * * * ****************** -,V 



SUBROUTINE OPTION( IOPT1 , IOPT2 , NUM ,RIRAY , R JRAY , RKRAY , RARAY , 

+ TARAY , AINRAY , APRAY , TIMRAY ) 

THIS SUBROUTINE GIVES THE USER A CHOICE OF OPERATIONS THAT CAN BE 
PERFORMED ON THE PROGRAM AND RETURNS THE USERS CHOICE WITH 
VARIABLES IOPT1 AND IOPT2 

DIMENSION RIRAY( 500) ,RJRAY(500) ,RKRAY(500) ,RARAY( 500) ,TARAY(500), 
+ A I NRAY ( 500) , APRAY( 500) ,TIMRAY( 500) 

CHARACTER* l,YORN 
IOPT1 = 0 

PROMPT USER FOR OPTION 

103 PRINT*, ’WHICH OF THE FOLLOWING OPTIONS WOULD YOU LIKE: ’ 

PRINT*,’ 1. -CALCULATE THE NEXT ORBIT USING THE SAME’ 

PRINT*,’ PARAMETERS' 

PRINT*,' 2. -CHANGE THE INITIAL PARAMETERS OF THE ORBIT' 

PRINT*,' 3. -CHANGE THE VELOCITY AT A POINT IN THE ORBIT' 

PRINT*,’ (THE UNPERTURBED ORBIT WILL BE USED)' 

PRINT*,' 4. -PLOT ANOTHER VIEW OF THE ORBIT(S)' 

PRINT*,' 5. -QUIT' 

PRINT*, 'ENTER 1, 2, 3, 4, OR 5:' 

READ* , IOPT2 

PRINT* , IOPT2 

CALL EXCMS( ’ CLRSCRN' ) 

IF ( IOPT2 .GT. 5) THEN 
GOTO 103 
END IF 

Prompt USER FOR TYPE OF ORBIT DESIRED 
105 IF ( IOPT2 . EQ. 1) THEN 

PRINT*, 'WHICH TYPE OF ORBIT WOULD YOU LIKE TO SEE,' 

PRINT*,' 1. -Unperturbed ORBITS’ 
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PRINT'-,' 2. -Perturbed ORBITS' 

PRINT* , ' ENTER 1 OR 2: ' 

READ* , IOPT1 

PRINT*, IOPT1 

CALL EXCMS('CLRSCRN') 

IF ( ( IOPT1 . NE. 1) .AND. (IOPT1 . NE. 2)) THEN 
PRINT*, ’ INVALID ENTRY! ' 

GOTO 105 

END IF 
END IF 

* PROMPT USER TO CLEAR PREVIOUS ORBITS 

107 IF ( ( I0PT2 .EQ. 1) .OR. (I0PT2 . EQ. 3)) THEN 

PRINT*, 'DO YOU WANT TO CLEAR THE PREVIOUS ORBITS?' 

PRINT*, 'ENTER "Y" OR "N" : ' 

READ*,YORN 

PRINT*, YORN 

CALL EXCMS( ' Clrscrn' ) 

IF (YORN .EQ. ' Y ' ) THEN 

RIRAY(l) = RIRAY(NUM) 

RJRAY(l) = RJRAY(NUM) 

RKRAY(l) = RKRAY(NUM) 

RARAY(l) = RARAY(NUM) 

TARAY(l) = TARAY(NUM) 

AINRAY(l) = AINRAY(NUM) 

A PRAY ( 1) = APRAY(NUM) 

TINRAY(l) = TIMRAY(NUM) 

NUM = 1 

ELSEIF (YORN . NE. ' N ’ ) THEN 
PRINT*, ' INVALID ENTRY! ! ' 

PRINT*, 'ALL INPUTS MUST BE CAPITOL LETTERS' 

GOTO 107 

ENDIF 
END IF 

* CHECK FOR INVALID OPTION 

IF ( ( IOPT2 . NE. 1). AND. (IOPT2 . NE. 2). AND. ( IOPT2 . NE. 3) .AND. 
+ ( IOPT2 .NE. 4). AND. (IOPT2 .NE. 5)) THEN 

PRINT*, ’ INVALID ENTRY!' 

GOTO 103 
ENDIF 
RETURN 
END 



***************************************************************** 
* COORDINATE TRANSFORMATIONS 

***************************************************************** 



SUBROUTINE PQWI JK( LAN , AP , INC , P ,Q , W , I , J ,K) 

THIS SUBROUTINE TRANSFORMS PQW COORDINATES TO IJK COORDINATES 

DOUBLE PRECISION INC , P ,Q, V, I , J ,K,R11 ,R12 ,R13 ,R21 ,R22 ,R23 , 

+ R31 ,R32 ,R33 , LAN, AP 

Rll = DCOS(LAN)*DCOS(AP) - DSIN( LAN)*DSIN( AP)*DCOS( INC) 

R12 = -DCOS(LAN)*DSIN(AP) - DSIN( LAN)*DCOS( AP)*DCOS( INC) 

R13 = DS IN( LAN )*DS IN( INC ) 
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R21 = DSI.N'C LAN)*DCOS( AP) + DCOS(LAN)*DSIN( AP)*DCOS( INC) 
R22 = -DSIN( LAN)*DSIN( AP) + DCOS( LAN)*DCOS( AP)*DCOS( INC) 
R23 = -DCOS( LAN)*DSIN( INC) 

R31 = DSIN( AP)*DSIN( INC) 

R32 = DCOS( AP)*DSIN( INC) 

R33 = DCOS(INC) 

I = R11*P + R 1 2*Q + R13*W 

J = R21*P + R22*Q + R23*W 

K = R3 1*P + R32*Q + R33*W 

RETURN 
END 



SUBROUTINE I JKPQW( LAN , AP , INC , I , J ,K , P , Q , W) 

THIS SUBROUTINE TRANSFORMS IJK COORDINATES TO PQV COORDINATES 

DOUBLE PRECISION INC , I , J ,K, P ,Q,W ,R1 1 ,R12 ,R13 ,R21 ,R22 ,R23 , 

+ R31 ,R32,R33 , LAN.AP 

Rll = DCOS( LAN) "DCOS( AP) - DSIN( LAN)*DSIN( AP)*DCOS( INC) 

R21 = -DCOS(LAN)*DSIN( AP) - DSIN(LAN)*DCOS( AP)*DCOS( INC) 

R31 = DSIN( LAN)*DSIN( INC) 

R12 = DSIN( LAN)*DCOS( AP) + DCOS(LAN)*DSIN( AP)*DCOS( INC) 

R22 = -DSIN( LAN) "DSIN( AP) + DCOS(LAN)*DCOS( AP)*DCOS( INC) 

R32 = -DCOS( LAN)*DSIN( INC) 

R13 = DSIN( AP)*DSIN( INC) 

R23 = DCOS( AP)*DSIN( INC) 

R33 = DCOS(INC) 

P = Rll* I + R 1 2* J + R13*K 

Q = R21*I + R22*J + R23*K 

W = R31*I + R32*J + R33*K 

RETURN 
END 






SUBROUTINE I JKRSW( LAN , AL , INC , I , J ,K , R , S , W) 

THIS SUBROUTINE CHANGES FROM IJK COORDINATES TO RSW COORDINATES 

DOUBLE PRECISION INC , I , J ,K,R, S , V ,R1 1 , R12 ,R13 ,R2 1 ,R22 ,R23 , 

+ R3 1 , R3 2 , R3 3 , LAN , AL 

Rll = DCOS( LAN)*DCOS( AL) - DSIN( LAN)*DCOS( INC)*DSIN( AL) 

R12 = DSIN( LAN)*DCOS( AL) + DSIN( AL)*DCOS( LAN)*DCOS( INC) 

R13 = DSIN( INC)*DSIN( AL) 

R21 = -DCOS(LAN)*DSIN( AL) -DSIN( LAN)*DCOS( INC)*DCOS(AL) 

R22 = -DSIN( LAN)*DSIN( AL) + DCOS( LAN)*DCOS( INC)*DCOS( AL) 

R23 = DSIN( INC)*DCOS( AL) 

R31 = DSIN( LAN)*DSIN( INC) 

R32 = -DCOS(LAN)*DSIN( INC) 

R33 = DCOS(INC) 

R = R1 1*1 + R12*J + R13*K 

S = R2 1*1 + R22*J + R23*K 

W = R3 1*1 + R32*J + R33*K 

RETURN 
END 
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SUBROUTINE RSWI JK( LAN, AL, INC ,R , S , W, I , J ,K) 

THIS SUBROUTINE CHANGES FROM RSW COORDINATES TO IJK COORDINATES 

DOUBLE PRECISION INC ,R , S ,W, I , J,K,R1 1 ,R12 ,R13 ,R2 1 , R22 ,R23 , 

+ R31 , R32 ,R33 ,LAN, AL 

RU = DCOS(LAN)*DCOS(AL) - DSIN( LAN)*DCOS( INC )*DSIN( AL) 

R21 = DSIN( LAN)*DCOS( AL) + DSIN( AL)*DCOS( LAN)*DCOS( INC ) 

R31 = DSIN( INC)*DSIN( AL) 

R12 = -DCOS( LAN) ,V DSIN( AL) -DSIN( LAN)*DCOS( INC )*DCOS( AL) 

R22 = -DSIN(LAN)*DSIN(AL) + DCOS( LAN)*DCOS( INC )*DCOS( AL) 

R32 = DSIN( INC)*DCOS( AL) 

R13 = DSIN(LAN)*DSIN( INC) 

R23 = -DCOS( LAN)*DSIN( INC ) 

R33 = DCOS(INC) 

I = R 1 1*R + R12*S + R13*W 

J = R21*R + R22"S + R23*W 

K = R31*R + R32 ,V S + R33*W 

RETURN 
END 
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SUBROUTINE PQWRSW(TA , P , Q , W , R , S , WN) 

THIS SUBROUTINE CHANGES FROM PQU COORDINATES TO RSW COORDINATES 

DOUBLE PRECISION P,Q,W,R, S ,WN,RU ,R12 ,R13 ,R21 ,R22 ,R23 , 

+ R31 ,R32 ,R33 ,TA 

Rll = DCOS(TA) 

R12 = DSIN(TA) 

R13 = 0.0 
R21 = -DSIN(TA) 

R22 = DCOS(TA) 

R23 = 0. 0 
R31 = 0. 0 
R32 = 0.0 
R33 = 1. 0 

R = R11*P + R 1 2*Q + R 1 3*W 
S = R2i*P + R22--Q + R23*W 
WN = R31--P + R32*Q +R33*W 
RETURN 
END 






SUBROUTINE RSWPQW(TA,R , S ,W ,P ,Q ,WN) 

THIS SUBROUTINE CHANGES FROM RSW COORDINATES TO PQW COORDINATES 

DOUBLE PRECISION R , S ,W,P ,Q,WN,R11 ,R12 ,R13 ,R21 ,R22 ,R23 , 

+ R31 ,R32 ,R33 ,TA 

Rll = DCOS(TA) 

R21 = DSIN(TA) 

R31 = 0. 0 
R12 = -DSIN(TA) 
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R22 = DCOS(TA) 

R32 = 0. 0 
R13 = 0. 0 
R23 = 0. 0 
R33 = 1. 0 

P = R11*R + R12*S + R13*W 
Q = R21*R + R22*S + R23*W 
WN = R31*R + R32*S +R33*W 
RETURN 
END 



* STORE ELEMENTS IN ARRAYS 
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SUBROUTINE STORE (RI,RJ,RK,R,TA,RI RAY , R JRAY , RKRAY , RARAY , TARAY , NUM , 
+ T AP ATNRAY APRAY TT TIMRAY"! 

THIS’ SUBROUTINE STORES THE* POSITION AND ELEMENTS IN ARRAYS IN 
SINGLE PRECISION FORM, FOR PLOTTING 

DOUBLE PRECISION RI ,RJ,RK,R ,TA , I , AP ,TT 

DIMENSION RIRAY(SOO) ,RJRAY(500) ,RKRAY(500) ,RARAY(500) ,TARAY(500) , 
+ AINRAYC500) ,APRAY(500) ,TIMRAY(500) 

R I RAY (NUM) = SNGL(RI) 

R JRAY (NUM) = SNGL(RJ) 

RKRAY (NUM) = SNGL(RK) 

RARAY (NUM) = SNGL(R) 

TARAY (NUM) = SNGL(TA) 

A INRAY (NUM) = SNGL(I) 

APRAY (NUM) = SNGL(AP) 

TIMRAY(NUM) = SNGL(TT) 

RETURN 

END 
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* INITIAL POSITION, VELOCITY 

Vr VrVrVr VrVrVrVrVrVrVr Vr VrVr VrVr VrVr VrVr Vr* Vr VrVcVr VrVr Vr Vr Vr Vr Vr VnV Vr Vr** * Vr Vr Vr VrVrVr VrVrVr VrVr VrVr VrVrVr Vr Vr VrVrVr VrVrVr Vr 



SUBROUTINE INPUTS( RI , RJ , RK , R , VI , VJ , VK , V , MU , QUIT , PI ) 

THIS SUBROUTINE GIVES THE USER A CHOICE TO EITHER ENTER THE 
INITIAL POSITION AND VELOCITY VECTOR OR TO LET THE PROGRAM 
CALCULATE THE INITIAL POSITION AND VELOCITY FROM USER PROMPTED 
INPUTS 

SUBROUTINES CALLED FROM THIS SUBROUTINE: 

INELTS = Prompts USER FOR ORBITAL ELEMENTS 

IPOS = PROMPTS USER FOR INITIAL POSITION (IJK) 

IVEL = PROMPTS USER FOR INITIAL Velocity (IJK) 

DOUBLE PRECISION RI ,RJ,RK,R,VI ,VJ,VK,V,MU,PI 
CHARACTER* 1, QUIT 

PROMPT USER FOR METHOD TO ENTER INPUTS 
195 PRINT*, ’IN WHICH MANNER WOULD YOU LIKE TO INPUT THE INITIAL’ 
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POSITION AND VELOCITY OF THE SATELLITE?' 

1: BY Inputting THE INITIAL POSITION AND VELOCITY* 
VECTORS IN THE PERIFOCAL COORDINATE SYSTEM (IJK)' 
2: BY LETTING THE SATELLITE BE PLACED ON THE "i"' 
AXIS OF THE (IJK) SYSTEM AT A DESIRED RADIUS OF' 
PERIGEE(RP) AND INPUTTING EITHER A DESIRED RADIUS 
OF APOGEE(RA) , A DESIRED ECCENTRIC ITY(E) , OR THE' 
DESIRED VELOCITY AT THAT RADIUS, AND A DESIRED' 
INCLINATION^ I). ’ 

3: QUIT' 

ENTER 1, 2 OR 3: ' 



PRINT-- , ’ 

FRINT* , ’ 

PRINT* , ’ 

PRINT*, ' 

PRINT* , ' 

PRINT*, ' 

PRINT*, ’ 

PRINT*, ' 

PRINT* , ' 

PRINT*, ' 

PRINT*, ' 

READ* , ICHC 

PRINT* , ICHC 

CALL EXCMS( 'CLRSCRN' ) 

• USER INPUTS POSITION AND VELOCITY VECTORS 
IF (ICHC .EQ. 1) THEN 

CALL IPOS(RI ,RJ,RK,R) 

CALL I VEL( VI , VJ , VK , V , R , MU) 

USER INPUTS ORBITAL ELEMENTS TO GET POSITION AND VELOCITY 
ELSEIF (ICHC .EQ. 2) THEN 

CALL INELTS(RI ,RJ,RK,R,VI ,VJ,VK,V,MU,PI) 

STOP PROGRAM 

ELSEIF (ICHC .EQ. 3) THEN 
QUIT = 'N' 

ELSE 

PRINT*, ' INVALID ENTRY! TRY AGAIN!' 

GOTO 195 
END IF 
RETURN 
END 

SUBROUTINE IPOS(RI ,RJ,RK,R) 

• THIS SUBROUTINE ASKS THE USER FOR THE INITIAL POSITION OF THE 

• Satellite IN GEOCENTRIC -EQUATORIAL COORDINATE SYSTEM 

DOUBLE PRECISION RI,RJ,RK,R 

CHARACTER*1 , CHOICE 
LOGICAL CORREC 
CORREC = .FALSE. 

' PROMPT USER FOR VELOCITY VECTOR 

180 IF(. NOT. CORREC) THEN 

CALL EXCMS( ’CLRSCRN' ) 

PRINT*, 'ENTER RADIUS VECTOR VALUES IN "KM"' 

PRINT*, 'RADIUS OF THE EARTH = 6400 KM' 

CORREC = .TRUE. 

PRINT*, 'ENTER RI : ' 

READ* , R I 

PRINT*, 'RI = ' ,RI , ' KM' 

PRINT*, 'ENTER RJ : ’ 
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READ' - ' , R J 

PRINT*, 'RJ = ' , R J , ' KM ' 

PRINT*, 'ENTER RK : ' 

READ* , RK 

PRINT*, 'RK = ’ ,RK, 'KM' 

CALCULATE TOTAL R 

R = DSQRT( (RI**2) + (RJ**2) + (RK**2)) 

PRINT*, 'R = ' ,R, 'KM' 

IF (R .LE. 6400.0) THEN 

PRINT*, 'RADIUS TO SMALL!! ENTER NEW VALUES!! ' 
GOTO 180 

END IF 

CHECK WITH USER THAT Values ARE CORRECT 
PRINT* , ' ARE THESE VALUES CORRECT?' 

PRINT*,' ENTER "Y" OR "N" 

READ*, CHOICE 
CHOICE = ' Y’ 

PR I NT*, CHOICE 
IF (CHOICE. EQ. ' Y' ) THEN 
CORREC = .TRUE. 

END IF 
GOTO 180 
END IF 
RETURN 
END 



SUBROUTINE IVEL( VI , VJ , VK , V , R , MU) 

THIS SUBROUTINE ASKS THE USER FOR THE INITIAL VELOCITY OF THE 
Satellite 

DOUBLE PRECISION VI , VJ , VK , V ,R , VCIR , VMAX ,MU 

CHARACTER* 1 , CHOICE 
LOGICAL CORREC 
CORREC = .FALSE. 

CALCULATE ESCAPE VELOCITY AND CIRCULAR VELOCITY AND PROMPT USER 
FOR VELOCITY VECTOR 
190 IF(. NOT. CORREC) THEN 

CALL EXCMS( 'CLRSCRN' ) 

VCIR = DSQRT(MU/R) 

VMAX = DSQRT((2. 0*MU)/R) 

PRINT*, 'CIRCULAR VELOCITY = VCIR KM/SEC ' 

PRINT*, 'MAXIMUM VELOCITY = ', VMAX, 'KM/ SEC* 

CORREC = .TRUE. 

PRINT*, 'ENTER VELOCITY VECTOR IN (KM/SEC)' 

PRINT*, 'ENTER VI : ' 

READ* , VI 

PRINT* , ' VI = ' , VI , ' KM/SEC ’ 

PRINT*, 'ENTER VJ : ' 

READ* , V J 
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PRINT*,' VJ = ' ,VJ, 'KM/SEC' 

PRINT*, 'ENTER VK : ' 

READ* , VK 

PRINT*,' VK = ' ,VK, 'KM/SEC' 

* CALCULATE TOTAL VELOCITY (V) 

V = DSQRT((VI**2) + ( VJ**2 ) + (VK**2)) 

PRINT*, 'V = ' ,V, 'KM/SEC' 

* CHECK WITH USER THAT VALUES ARE CORRECTS 
PRINT-*, 'ARE THESE VALUES CORRECT?' 

PRINT*,' ENTER "Y" OR "N" :' 

RE AD*, CHOICE 
CHOICE = ' Y’ 

PRINT*, CHOICE 
IF (CHOICE. EQ. ' Y' ) THEN 
CORREC = .TRUE. 

END IF 

IF (V . GE. VMAX) THEN 

PRINT*, 'VELOCITY IS GREATER THAN THE ESCAPE VELOCITY!! ' 
PRINT*, 'RE-ENTER VELOCITY! ! ! ' 

CORREC = .FALSE. 

ENDIF 
GOTO 190 
ENDIF 
RETURN 
END 
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SUBROUTINE INELTS(RI , RJ , RK, R, VI, VJ, VK, V,MU,PI) 

SATELLITE PLACED ON ' I ' AXIS AND USER SUPPLY ORBITAL ELEMENTS TO 
GET INITIAL POSITION AND VELOCITY 

DOUBLE PRECISION RI , RJ ,RK , R , VI , VJ , VK , V ,MU , I , ENR ,A ,E ,RP ,RA ,PI , VMAX 
CHARACTER* 1, CHOICE 

* PROMPT USER FOR PERIGEE RADIUS 

198 PRINT*, ' ENTER RADIUS OF PERIGEE(RP) IN (KM), FOR EXAMPLE:' 

PRINT*, 'LOW EARTH ORBIT (LEO), RP = 6600.0 KM' 

PRINT*, 'GEOS YNCROUNOUS ORBIT, RP = 42241.1 KM' 

PRINT*, 'ENTER RP: ' 

PRINT* , ' "RP" MUST BE > 6400KM' 

READ* , RP 
PRINT*, RP 

* CHECK FOR VALID RADIUS 
IF (RP . LT. 6400.0) THEN 

PRINT*, 'YOUR "RP" IS TO SMALL!!' 

GOTO 198 
ENDIF 

* PROMPT USER FOR TYPE OF INPUT 

PRINT**, 'DO YOU WANT TO ENTER THE ECCENTRICITY (E), ' 

PRINT*, 'RADIUS OF APOGEE (RA), OR VELOCITY (V)?' 

PRINT*, 'ENTER "E", "R", OR "V": ' 
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READ*, CHOICE 

PR I NT*, CHOICE 

CALL EXCMS( 'CLRSCRN' ) 

USER ENTERS Eccentricity AND SEMI -MAJOR AXIS, ENERGY AND VELOCITY 
IS CALCULATED IN THAT ORDER 
IF (CHOICE .EQ. ’ E ' ) THEN 

PRINT** , ' ENTER ECCENTRICITY (E): ' 

PRINT**, '0. 0 <= E < 1. O' 

READ** , E 
PRINT* ,E 

CHECK FOR VALID ECCENTRICITY 
IF ( ( E .LT. 0.0) .OR. (E . GE. 1.0)) THEN 
PRINT*,' INVALID "E"' 

GOTO 198 
END IF 

A = RP/(1-E) 

ENR = -MU/( 2. 0*A) 

V = DSQRT( 2*(ENR+(MU/RP) ) ) 

USER INPUTS RADIUS OF APOGEE AND ECCENTRICITY' IS CALCULATED 
THEN SEMI -MAJOR AXIS, ENERGY AND THEN VELOCITY. 

ELSEIF (CHOICE . EQ. ’ R* ) THEN 

PRINT*, 'ENTER RADIUS OF APOGEE (RA) IN KM: ' 

PRINT*, '"RA" MUST BE >="RP", "RP" = ' ,RP 
READ* , RA 
PR I NT*, R A 

CHECK FOR VALID RADIUS OF APOGEE 
IF (RA .LT. RP) THEN 

PRINT*, 'YOUR "RA" IS TO SMALL!!' 

GOTO 198 
END IF 

E = ( RA-RP) / ( RA+RP) 

A = RP/( 1-E) 

ENR = -MU/ (2. 0*A) 

V = DSQRT( 2*(ENR+(MU/RP) ) ) 

USER INPUTS MAGNITUDE OF VELOCITY, PROGRAM PROVIDES CIRCULAR 
AND ESCAPE VELOCITY FOR COMPARISON AND TO CHECK FOR VALID 
INPUTS 

ELSEIF (CHOICE . EQ. ' V ' ) THEN 

PRINT*, ’ENTER VELOCITY IN KM/SEC: ’ 

PRINT*, ’THE MINIMUM VELOCITY ALLOWED IS FOR A CIRCULAR ORBIT’ 
VC IRC = SQRT( SNGL(MU/RP) ) 

PRINT*, ’ ORBIT. V( Circular) = ’.VCIRC,’ KM/S’ 

VMAX = DSQRT(2*(MU/RP)) 

PRINT**, ’THE MAXIMUM VELOCITY < ’,VMAX,’ KM/S' 

READ* ,V 
PRINT* ,V 

IF (V .LT. VCIRC) THEN 

PRINT*, ’VELOCITY’ TO SMALL! ’ 

GOTO 198 

END IF 

IF (V .GE. VMAX) THEN 
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PRINT* , ' VELOCITY TO GREAT!! ' 

GOTO 198 

ENDIF 

ELSE 

PRINT* INVALID ENTRY! TRY AGAIN' 

GOTO 198 
ENDIF 

INCLINATION NEEDED TO GIVE Velocity A Direction 
PRINT*, 'ENTER INCLINATION (I) IN DEGREES:' 
READ*, I 
PRINT*, I 

I = (PI/180. 0)*I 
VK = V*DSIN( I) 

VJ = V*DCOS( I) 

VI = 0. 0 

RADIUS VECTOR SET 

RI = RP 

RJ = 0. 0 

RK = 0. 0 

R = RP 

RETURN 

END 
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* CALCULATE THE ORBITAL ELEMENTS 



SUBROUTINE CALCEL( RI , RJ , RK , R , VI , VJ , VK , V , E I , E J , EK , E , A , I , LAN , 

+ LP ,TA , PER ,EA ,MA , AP , AL,TF , P , PI ,MU ,MM ,N , H,HI ,HJ) 

THIS SUBROUTINE CALLS THE INDIVIDUAL SUBROUTINES TO CALCULATE THE 
ORBITAL ELEMENTS 

THIS SUBROUTINE CALLS THE FOLLOWING SUBROUTINES( RETURNED VALUES) 
ENERGY = ENERGY PER MASS (ENR) 

ANGMOM = ANGULAR MOMENTUM (H,HI,HJ,HK) 

NODE = NODE VECTOR (N,NI ,NJ,NK) 

LATREC = SEMI -LATUS RECTUS (P) 

ECC = ECCENTRICITY (E,EI,EJ,EK) 

SMAXIS = SEMI -MAJOR AXIS (A) 

INCL = INCLINATION (I) 

ASNODE = LONGITUDE OF ASCENDING NODE (LAN) 

ARP = ARGUMENT OF PERIGEE (AP) 

IJKPQW = 'IJK' SYSTEM TO 'PQW' SYSTEM 
TANOM = TRUE ANOMALY (TA) 

ARLAT = ARGUMENT OF LATITUDE (AL) 

LONPER = LONGITUDE OF Perigee (L?) 

TLON = TRUE LONGITUDE (TL) 

PERIOD = PERIOD (PER) 

ECCAN = ECCENTRIC ANOMALY (EA) 

MEANMO = MEAN MOTION (MM) 

MSANAN = MEAN ANOMALY (MA) 

TFLGHT = TIME OF FLIGHT (TF) 

DOUBLE PRECISION RI ,RJ,RK,R, VI , VJ , VK, V ,EI ,EJ,EK, E , A, I , LAN , AL, 
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+ LP ,TA , PER ,EA ,MA . AP ,TF ,HI ,HJ ,HK ,H ,NI ,NJ, NK ,N , P , PI , MU , MM, ENR , 

+ TL , R? , RQ , RW , NP , NQ , NV 



CALL ENERGY ( V , R , MU , EN'R ) 

CALL ANGMOM( RI , RJ , RK , VI , V J , VK , HI , HJ , HK , H ) 

CALL NODE ( H I , H J , N I , N J , NK , N ) 

CALL LATREC(H,P,MU) 

CALL ECC(RI ,RJ ,RK,R , VI , VJ, VK, V ,EI ,EJ ,EK ,E ,MU) 

CALL SMAXI S ( MU , ENR , A) 

CALL INCL( HK, H , I , PI ) 

SPECIAL CASE IF INCLINATION = 0. 0 
IF (I.NE. 0.0) THEN 

CALL ASNODE ( NI , N , LAN , N J , P I ) 

CALL ARP( NI ,NJ,N ,EI ,EJ,EK,E , AP , PI ,NP ,NQ , LAN) 

ELSE 

LAN = 0. 0 
A? = 0. 0 
END IF 

COORDINATE TRANSFORMATION OF ' R’ AND V VECTORS 
CALL IJKPQW( LAN , AP , I , RI , R J , RK , RP , RQ , RW ) 

CALL I JKPQW( LAN , AP , I , N I , NJ , NK , NP , NQ , NV ) 

CALL TANOM(EI ,EJ,EK,E ,RI ,RJ,RK,RP,RQ,RV,R,VI ,VJ, VK,TA,PI) 

SPECIAL CASE FOR Inclination =0.0 
IF (I . NE. 0. 0) THEN 

CALL ARLAT(NI ,NJ,NK,N,RI ,RJ,RK,R,AL,PI ,TA,AP) 

ELSE 

AL = TA 
END IF 

CALL LONPERC LAN , AP , LP) 

CALL TLON ( LAN , AP , TA , TL ) 

CALL PERIOD( A , PER , PI ,MU) 

CALL ECCAN( E ,TA ,EA , PI ) 

CALL ME ANMO ( A , MM , MU ) 

CALL ME ANAN ( E A , E , MA ) 

CALL TFLGHT( MM , MA , TF ) 

RETURN 

END 
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SUBROUTINE ENERGY ( V , R , MU , ENR) 

THIS SUBROUTINE CALCULATES THE ENERGY OF THE ORBIT 

DOUBLE PRECISION V,R,MU,ENR 

ENR = ( ( V**2 ) / 2 ) - (MU/R) 

RETURN 

END 
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SUBROUTINE ANGMOM( RI , RJ , RK , VI , VJ , VK , HI , HJ , HK , H) 
THIS SUBROUTINE CALCULATES THE ANGULAR MOMENTUM 
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DOUBLE PRECISION' RI ,RJ,RK, VI , VJ, VK.HI ,HJ, HK ,H 

HI = (RJ * VK) - (RK * VJ) 

HJ = (RK * VI) - (RI * VK) 

HK = (RI * VJ) - (RJ * VI) 

H = DSQRT((HI**2) + (HJ**2) + (HK**2)) 

RETURN 

END 
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SUBROUTINE NODE (HI ,HJ,NI ,NJ,NK,N) 

THIS SUBROUTINE CALCULATES THE NODE VECTOR 

DOUBLE PRECISION HI ,HJ,NI ,NJ,NK,N 

NI = -HJ 
NJ = HI 
NK = 0. 0 

N = DSQRT((NI**2) + (NJ**2)) 

RETURN 

END 
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SUBROUTINE LATREC(H,P,MU) 

THIS SUBROUTINE CALCULATES THE SEMI-LATUS RECTUM 

DOUBLE PRECISION H,P,MU 

P = ( H**2 ) /MU 

RETURN 

END 
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SUBROUTINE ECC(RI ,RJ,RK,R, VI ,VJ,VK,V,EI ,EJ,EK,E,MU) 

* THIS SUBROUTINE CALCULATES THE ECCENTRICITY 

DOUBLE PRECISION RI ,RJ,RK , R, VI , VJ, VK, V,EI ,EJ,EK,E , MU, DOT 

* CALCULATE DOT PRODUCT OF ' R* AND 'V' VECTORS 
DOT = (RI*VI) + (RJ*VJ) + (RK*VK) 

El = ( 1. OD+OO/MU) * ( ( ( V ,v *2) - (MU/R)) * RI - (DOT)*VI) 

EJ = (1. OD+OO/MU) (((V**2) - (MU/R)) * RJ - (DOT)*VJ) 

EK = (1. OD+OO/MU) * ( ( ( V**2 ) - (MU/R)) * RK - (DOT)*VK) 

E = DSQRT( (EI**2) + (EJ**2) + (EK**2)) 

RETURN 

END 






SUBROUT I NE SMAX I S ( MU , ENR , A ) 

* THIS SUBROUTINE Calculates THE SEMI -MAJOR AXIS 
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DOUBLE PRECISION MU , ENR , A 



A = -MU/( 2 -ENR) 

RETURN 

END 



*■ 
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SUBROUTINE INCL(HK ,H , I , PI ) 

THIS SUBROUTINE CALCULATES THE INCLINATION 
'I' ALWAYS LESS THAN 180 DEGREES 

DOUBLE PRECISION HK,H,I,PI 

I = DACOS( HK/H) 

RETURN 

END 
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SUBROUTINE ASNODE(NI , N , LAN , NJ , PI ) 

THIS SUBROUTINE CALCULATES THE LONGITUDE OF THE ASCENDING NODE 
IF 'NJ' > 0 THEN 'LAN' < 180 DEGREES 

DOUBLE PRECISION NI ,N,LAN,NJ,PI 



LAN = DATAN2(NJ ,NI ) 

IF (LAN . LT. 0.0) THEN 
LAN = ( 2 "PI ) + LAN 
END IF 
RETURN 
END 






JU 






SUBROUTINE ARP( NI , NJ , N , El , E J , EK , E , AP , PI , NP , NQ , LAN) 

THIS SUBROUTINE CALCULATES THE ARGUMENT OF Perigee 

IF ’EK' GREATER THAN 0 THEN 'AP' < 180 

VARIABLE TEMP USED AS A Temporary VALUE FOR ARCTAN 

DOUB LE PRECISION NI , N J , N , E I , E J , EK , E , AP , P I , NQ , NP , TEMP , LAN 

IF ((El . EQ. 0.0) .AND. (EJ . EQ. 0.0)) THEN 
AP = 0. 0 
ELSE 

TEMP = DATAN2(EJ,EI) 

IF (TEMP . GT. LAN) THEN 
AP = TEMP - LAN 

ELSE 

AP = (2*PI) - (LAN - TEMP) 

ENDIF 

IF ( AP .LT. 0.0) THEN 
AP = (2*PI) + AP 
ENDIF 

IF (AP . GT. ( 2*PI) ) THEN 
AP = A? - (2*PI) 
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END IF 
END IF 
RETURN 
END 



s{*t*%*\** T ^ f 






SUBROUTINE TANOM( E I , E J , EX , E , RI , RJ, RX , R? , RO , RV 
+ TA , ?I ) 

THIS SUBROUTINE CALCULATES THE TRUE Anomaly 
IF (R DOT V) > 0 THEN TA < 180 DEGREES 



V I , V J , VK , 



DOUBLE PRECISION 
+ RP.RQ.RV 



DOT, El , EJ,EX,E ,RI ,RJ ,RX,R, VI , VJ , VX,' 



TA = DATAN2(RQ,RP) 

IF (TA . LT. 0. 0 ) THEN 
TA = (2 * PI) -r TA 
END IF 
RETURN 
END 



SUBROUTINE ARLAT ( NI , NJ , NX , N , RI , RJ . RX , R , AL , PI , TA 
THIS SUBROUTINE CALCULATES THE ARGUMENT OF LATI 
IF (RX > 0) THEN AL < 150 DEGREES 




DOUBLE PRECISION NI ,NJ,NX,N,R! 



ST iS 

. ~ y J 



AL = TA 



RETURN 



END 
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SUBROUTINE LONPER( LAN , A? , LP) 

THIS SUBROUTINE CALCULATES THE LONGITUDE OF PERIGEE 
DOUBLE PRECISION LAN , A? , L? 



L? = LAN + A? 

RETURN 

END 



SUBROUTINE TLON( LAN , AP , TA , TL) 

THIS SUBROUTINE CALCULATES THE TRUE LONGITUDE AT EPOCH 

DOUBLE PRECISION LAN , A? , TA , TL 

TL = AP + LAN + TA 

RETURN 

END 
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SUBROUTINE PERIOD(A,PER,PI ,MU) 

* THIS SUBROUTINE CALCULATES THE PERIOD 

DOUBLE PRECISION A, PER, PI, MU 

PER = 2. 0D+00*(PI)*DSQRT((A**3)/MU) 

RETURN 

END 
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SUBROUTINE ECCAN(E,TA,EA,PI) 

* THIS SUBROUTINE CALCULATES THE ECCENTRIC Anomaly 

DOUBLE PRECISION E,TA,EA,PI 

EA = DACOS( (E + DCOS(TA) )/( 1. OD+OO + E*DCOS(TA))) 
IF (TA . GT. PI) THEN 
EA = (2*PI) - EA 
END IF 
RETURN 
END 
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SUBROUTINE MEANMO( A , MM , MU) 

THIS SUBROUTINE CALCULATES THE MEAN MOTION 

DOUBLE PRECISION A, MM, MU 

MM = DSQRT(MU/(A**3)) 

RETURN- 

END 
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SUBROUTINE MEANAN( EA , E , MA) 

* THIS SUBROUTINE CALCULATES THE MEAN Anomaly 

DOUBLE PRECISION EA,E ,MA 

MA = EA - E*DSIN(EA) 

RETURN 

END 
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SUBROUTINE TFLGHT ( MM , MA , TF ) 

THIS SUBROUTINE CALCULATES THE TIME OF FLIGHT 

DOUBLE PRECISION MM , MA , TF 

TF = (1/MM)*MA 
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end 



CALCULATE UNPERTURBED ORBIT 



SUBROUTINE UNPRET( DT , PER , AL , LAN , A? , I , R I , R J , RI\ , R , 

+ VI VJ V'< V M U - * - V Ti p w \<x ~ x 

4- TF ,T, NUM , R I RAY , RJRAY , RKRAY , RARAY , TARAY , A I NRAY , ARRAY , T I MHA' 

+ i .) 

THIS SUBROUTINE CALCULATE THE UNPERTURBED CRB IT 



THIS SUBROUTINE CALLS THE FOLLOWING SUBROUTINES: 
NEWELT = CALCULATE NEW ELEMENTS AFTER TIME STEP 
NEWPOS = CALCULATE NEW POSITION AFTER TIME STEP 
NEWVEL = CALCULATE NEW VELOCITY AFTER TIME STEP 
STORE = STORES POSITION IN ARRAYS 



DOUBLE 

+ 



PRECISION 



T 

1 j -* - i 




> "* L y — -"‘N , 

PI.H.A.E 








, R J , RI\ , R , V I , V J , UK 
, MM , MA , EA , TF , TT 



DIMENSION 



RARAY( 500 ) , TARAY( 500 ) , RIRAY( 500 ) , RJRAY( 500 ) , 
RKRAY(SOO) ,AINRAY(500) ,APRAY(500) ,TIMRAY(50C) 



SET TRUE ANOMALY TO NEGATIVE SO LOO? CAN BE EXECUTED 
IF (TA . GT. 6.21) THEN 
TA = TA - ( 2--PI) 

END IF 



CONTINUE AROUND CEBIT TILL CLOSE TO PERIGEE 
230 IF ((TA . LE. 6.21) .AND. (T . LE. PER) ) THEN 

Increment TRUE TIME 
TT = TT - DT 

CALL NEWELT ( MM , MA , E , EA , TA , TF , DT , PI . PER) 

CALL NPOS ( RI , RJ , RX , R , LAN , AP , I , TA , A , E ) 

CALL NVEL(E , P , TA , LAN , A? , I , VI , VJ , VK , V , MU) 

INCREMENT STEP COUNTER AND STORE VALUES 
NUM = NUM + 1 

CALL STORE( RI , RJ , RK . R , TA . RIRAY , RJRAY . RKRAY , 

+ RARAY , TARAY , NUM , I , AP , A I NRAY , APRAY , 

+ TT.TIMRAY) 



INCREMENT TIME STEP COUNTER 
T= T + DT 
GOTO 230 
END IF 
RETURN 
END 



* CALCULATE THE UNPERTURBED NEW ELEMENTS 
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SUBROUTINE NEWELT( MM , MA , E , EA , TA ,TF , DT , PI , PER ) 

THIS SUBROUTINE CALCULATES THE Unperturbed NEW ELEMENTS 

THIS SUBROUTINE CALLS THE FOLLOWING SUBROUTINES: 

NEA = NEW ECCENTRIC ANOMALY 
NT A = NEW TRUE ANOMALY 

DOUBLE PRECISION MM,MA,E,EA,TA,TF,DT,PI ,PER 

Increment TIME OF FLIGHT AND CHECK IF TF GREATER THAN PERIOD 
TF = TF + DT 
IF (TF . GT. PER) THEN 
TF = TF - PER 
END IF 

CALCULATE MEAN ANOMALY AND USE TO FIND ECCENTRIC Anomaly THEN NEW 
TRUE ANOMALY 
MA = MM*(TF) 

CALL NEA(MA,E ,EA) 

CALL NTA(EA,E,TA,PI) 

RETURN 

END 



CALCULATE PERTURBED ORBIT 
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SUBROUTINE PRETURC DT , PER , AL , LAN , AP , I , R I , RJ , RK , R , 

+ VI , VJ , VK , V ,FR ,FS ,FW ,MU, PI ,H , A ,E ,N,TA , P ,MM ,MA,EA , 

+ TF , T , NUM , RIRAY , RJRAY , RKRAY , RARAY , TAR AY , A INRAY , APR AY , T I MR AY , 

+ TT , TFEA , TFSU , TFMO , TFDRA , TDI , TDA , TDE , TDMM , TDMA , TDLAN , TDH , TDAP) 

THIS SUBROUTINE CALCULATES THE PERTURBED ORBIT. 

THIS SUBROUTINE CALLS THE FOLLOWING SUBROUTINES: 

TFORCE = CALCULATE THE TOTAL PERTURBING FORCE ON THE SATELLITE 

PNEWEL = CALCULATE THE Perturbed NEW ELEMENTS 

NPOS = NEW POSITION AFTER TIME STEP 

NVEL = NEW VELOCITY AFTER TIME STEP 

PERIOD = PERIOD OF PERTURBED ORBIT 

STORE = STORE POSITION AND ELEMENTS IN ARRAYS FOR PLOTTING 

DOUBLE PRECISION T,DT, PER, AL, LAN, AP , I ,RI ,RJ,RK, R,VI , VJ, VK,V, 

+ FR,FS ,FW,MU,PI ,H,A,E ,N,TA,P,MM,MA,EA,TF,TT, 

+ DI,DA,DE,DMM,DMA,DLAN,DH,DAP,EI,EJ,EK,HI,HJ,LP,M, 

+ DVR,DVS ,DVW,DVI ,DVJ,DVK 

DIMENSION RARAY(500) ,TARAY(500) ,RIRAY(500) ,RJRAY(500) , 

+ RKRAY(500) ,AINRAY(500) ,APRAY(500) ,TIMRAY(500) 

SET MEAN RADIUS OF EARTH 
RE = 6400.0 

DT = PER/50 
T = DT 

IF (TA .GT. 6.21) THEN 
TA = TA - (2 •■-PI) 
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END IF 

IF (TF . GE. PER) THEN 
TF = TF - FER 
EN'DIF 

CONTINUE Around ORBIT FOR ONE PERIOD 
IF ((TF . LT. PER) .AND. (T . LT. PER)) THEN 

INCREMENT TRUE TIME 
TT = TT + DT 

CALL TFORCE( AL , LAN , AP , I ,RI ,RJ,RK ,R , VI , VJ, VK, V, 

+ TT,FR,FS ,FW,MU,PI , 

+ FEA , FSU , FMO , FDRA , FOR , 

+ El ,EJ,EK,E ,A,T,LP,TA, PER,EA,MA,TF,P, 

+ MM,N,H,HI ,HJ,DT) 

CALL PNEWEL(FR,FS,FV,H,R,A,E,N,TA,DT, I ,LAN,AL, 

AP,P,MM,MA,EA,TF,T,MU,PI, 

DI , DA , DE , DMM , DMA , DLAN , DH , DAP ) 

CALL NP0S(RI , RJ,RK,R,LAN,AP,I , TA,A,E) 

CALL NVEL( E , P,TA , LAN, AP , I , VI , VJ, VK, V ,MU) 

CALCULATE NEW PERIOD AND RESET TIME STEP AND TIME COUNTER 
IF NOT AT END OF ORBIT 
IF (T .LT. (PER-DT)) THEN 
CALL PERIOD( A ,PER , PI ,MU) 

DT = PER/50 
T = TF 
END IF 

INCREMENT STEP COUNTER 
NUM = NUM + 1 

CALL STORE ( R I , R J , RK , R , TA , R I RAY , RJRAY , RKRAY , 

+ RARAY , TARAY , NUM , I , AP , AINRAY , APRAY , 

+ TT,TIMRAY) 

TOTAL ELEMENT CHANGES 
TDI = TDI + SNGL( ABS(DI) ) 

TDA = TDA + SNGL(ABS(DA) ) 

TDE = TDE + SNGL( ABS(DE) ) 

TDMM = TDMM + SNGL( ABS(DMM) ) 

TDMA = TDMA + SNGL( ABS(DMA) ) 

TDLAN = TDLAN + SNGL(ABS(DLAN) ) 

TDH = TDH + SNGL(ABS(DH)) 

TDAP = TDAP + SNGL(ABS(DAP) ) 

TFEA = TFEA + FEA 
TFSU = TFSU + FSU 
TFMO = TFMO + FMO 
TFDRA = TFDRA + FDRA 

CHECK FOR IMPACT 
IF (R . LE. RE) THEN 

PRINT* , 1 SATELLITE WILL IMPACT THE EARTH!! ' 

T = PER 

END IF 

INCREMENT TIME COUNTER 
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T = T + DT 
GOTO 240 
END IF 
RETURN 
END 



******* 



******************************************************* 



* CALCULATE THE PERTURBING FORCES 



* * * * * * * * * * * * * * * * * * * * * * ' 



****************************************** 



SUBROUTINE TFORCE( AL, LAN , AP , I ,RI ,RJ ,RK, R , VI , VJ , VK, V,TT, 

+ FR , FS , FW , MU , P I , FEA , FSU , FMO , FDRA , FOR , 

+ El ,EJ ,EK ,E , A ,T,LP ,TA ,PER , EA,MA ,TF , P , 

+ MM,N,H,HI,HJ,DT) 

THIS SUBROUTINE SUMS ALL THE PERTURBING FORCES FOR THE TOTAL 
PERTURBING FORCE. 

THE FOLLOWING SUBROUTINES WERE CALLED: 

OBERT = OBLATENESS OF THE EARTH 
FSUN = GRAVITATIONAL Attraction OF THE SUN 

FMOON = GRAVITATIONAL Attraction OF THE MOON 

FDRAG = DRAG FORCES 

DOUB LE PREC I SION FER , FES , FEW , FSR , FSS , FSW , FMR , FMS , FMW , MU , PI , 

+ FDR,FDS , FDW,FR,FS ,FW ,RI , RJ,RK,R, AL, I ,TT, LAN , AP , VI , VJ, VK, V, 

+ E I , E J , EK , E , A , T , LP , TA , PER , EA , MA , TF , P , 

+ MM,N,H,HI ,HJ,DT 

CALL OBE ART( RI , R J , RK , R , AL , I , FER , FES , FEW , MU ) 

CALL FSUN(TT,RI ,RJ,RK,R,FSR,FSS , FSW, PI) 

CALL FMOON ( TT , RI , RJ , RK , R , FMR , FMS , FMW , P I ) 

CALL FDRAG( RI , RJ ,RK ,R , VI , VJ, VK , V, LAN , AP , I , FDR ,FDS , FDW, 

+ El ,EJ ,EK , E , A,T, LP ,TA , PER,EA ,MA, AL,TF , P , PI ,MU, 

+ MM,N,H,HI ,HJ,DT) 



SUM VECTOR FORCES 
FR = FER + FSR + FMR + FDR 

FS = FES + FSS + FMS + FDS 

FW = FEW + FSW + FMW + FDW 



CALCULATE TOTAL FORCE FROM EACH, AND TOTAL OF ALL 
FEA = SNGL( SQRT( ( FER**2 ) +( FES**2 ) +( FEW**2 ) ) ) 

FSU = SNGL(SQRT((FSR**2)+(FSS**2)+(FSW**2))) 

FMO = SNGL( SQRT( ( FMR**2 ) +( FMS**2 ) +( FMW**2 ) ) ) 

FDRA = SNGL( SQRT( ( FDR**2 )+( FDS**2 )+( FDW**2 ) ) ) 

FOR = SNGL( SQRT( (FR**2)+(FS**2)+(FW**2) ) ) 

RETURN 

END 



****** 



’ * * ** * * * * * ** * * * * ** ’ 



**************************************** 



SUBROUTINE OBEART( RI , RJ , RK , R , AL , I , FER , FES , FEW , MU) 

THIS SUBROUTINE CALCULATES THE PERTURBING FORCE DUE TO THE 
OBLIQUENESS OF THE EARTH. 
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DOUBLE PRECISION J2 ,RE ,FER , FES , FEW , RI ,RJ ,RK, R , AL, I , MU ,M 

J2 = 1. 032364D-03 
RE = 6. 37820+03 

FER = ((-3. OD+O 0*MU* J 2* ( RE** 2 ) ) / ( 2 . 0D+00*(R**4) ))* 

+ (l.OD+OO - (3. OD+OO* ( (DSIN( I ) )**2) * ( (DSIN( AL) )**2) ) ) 

FES = ( ( - 3. OD+OO *MU* J 2* ( RE**2) ) / ( R**4 ) )* 

+ ( ( (DSIN( I) )**2)* (DSIN(AL))*(DCOS(AL))) 

FEW = ((-3. 0D+00*MU*J2*(RE**2))/(R**4))* 

+ (DSIN(I) *DCOS( I ) *DS I N( AL) ) 

RETURN 

END 

**************************************************************** 

SUBROUTINE FSUN(TT,RI ,RJ,RK,R,FSR,FSS,FSW,PI) 

* THIS SUBROUTINE CALCULATES THE PERTURBING FORCE DUE TO THE SUN 

* THE FOLLOWING SUBROUTINES ARE CALLED: 

SUNPOS = SUNS POSITION ORBITING AROUND EARTH 

* HEVBOD = PERTURBING FORCE FROM A Heavenly BODY 

DOUBLE PRECISION FSR.FSS ,FSW,PI , 

+ RSI , RSJ, RSK , S LAN, SI , SAL, SMU ,TT,RI , RJ ,RK, R ,RS 

* SUNS GRAVITATIONAL PARAMETER 
SMU = 1. 3271544D X 11 

CALL SUNPOS (TT , RSI , RSJ , RSK , RS , SLAN , S I , SAL , PI ) 



RETURN 

END 


















SUBROUTINE FMOON(TT , RI , R J , RK , R , FMR , FMS , FMW , PI ) 

THIS SUBROUTINE CALCULATES THE PERTURBING FORCE DUE TO The MOON 

THE FOLLOWING SUBROUTINE ARE CALLED: 

MONPOS = MOONS POSITION ORBITING AROUND THE EARTH 
HEVBOD = PERTURBING FORCE FROM A HEAVENLY BODY 

DOUBLE PRECISION FMR , FMS , FMW , RMI , RMJ , RMK , MLAN , MI , MAL , MMU , 

+ TT,RI,RJ,RK,R,RM,PI 

MOONS GRAVITATIONAL PARAMETER 
MMU = 4. 90287D+03 



CALL MONPOS (TT , RMI , RMJ , RMK , RM , MLAN , MI , MAL , PI ) 
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CALL HEVBOD( RI , RJ , RK , R , RMI , RMJ , RMK , RM , MLAN , MAL , MI , MMU , FMR , FMS , FMW) ORB 14430 



RETURN 

END 

SUBROUTINE HEVBODC RI , R J , RK , R , RP I , RPJ , RPK , RP , LAN , AL , INC , MUP , 
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+ FHR,FHS,FHW) 

THIS SUBROUTINE CALCULATES THE PERTURBING FORCE DUE TO A 
HEAVENLY BODY. 

THE FOLLOWING SUBROUTINE WAS CALLED: 

IJKRSW = 'IJK' SYSTEM TO THE 'RSW' SYSTEM 

DOUBLE PRECISION DOT,FHI ,FHJ,FHK,RI ,RJ,RK,R,RPI ,RPJ,RPK,RP, 
+ LAN, AL, INC ,MUP , I , J,K, IP , JP ,KP ,M,FHR,FHS ,FHW 



I = RI/R 
J = RJ/R 
K = RK/R 
IP = RPI/RP 
JP = RPJ/RP 
KP = RPK/RP 

CALCULATE DOT PRODUCT OF UNIT VECTORS 
DOT = (( I*IP )+( J*JP )+( K ,V KP )) 

CALCULATE FORCES IN THE ’IJK’ SYSTEM 
FHI = ( MUP/ ( RP**2) )*( R/RP)*( 3. 0D+00*D0T*( IP) -( I) ) 
FHJ = (MUP/(RP**2))*(R/RP)*(3. 0D+00*DOT*( JP) -( J) ) 
FHK = ( MUP/( RP**2) )*(R/RP)*( 3. OD+O0*DOT*( KP) -( K) ) 

Transform FORCES TO THE RSW SYSTEM 

CALL I JKRSW( LAN , AL , INC , FHI , FHJ , FHK , FHR , FHS , FHW) 

RETURN- 

END 



J- .'.J. J. J*. .1. J* J* .1. J. 



SUBROUTINE SUNPOS(TT , RSI , RSJ , RSK , RS , SLAN , SI , SAL, PI ) 
THIS SUBROUTINE CALCULATES THE SUNS POSITION 

VARIABLES USED TO DESCRIBE THE SUNS ORBIT: 

SI = SUNS INCLINATION 

SLAN= SUNS Longitude OF ASCENDING NODE 

SAP = SUNS ARGUMENT OF PERIGEE 

RS = SUNS ORBITAL RADIUS 

STA = SUNS TRUE ANOMALY 

SAL = SUNS ARGUMENT OF LONGITUDE 

DOUBLE PRECISION SLAN, SI , SAL,RS , STA, SAP ,TT, RSI ,RSK, 
+ RSJ.RSP ,RSQ,RSW,PI 

SI = 4. 09279709D-01 
SLAN = 0. 0D+00 
SAP = 0.0D+00 
RS = 1. 4959965D+08 

STA = ( ( 2. 0*PI )/( 365. 0 * 86400.0) * IT) 

SAL = STA + SAP 

CALCULATE SUNS POSITION IN 'PQW' SYSTEM 
RS? = RS-DCOS(STA) • 
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